Mercurial > repos > public > sbplib

comparison +scheme/AdvectionRV1D.m @ 1010:f753bada1a46 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add 1D scheme for advection with RV
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 07 Nov 2018 10:03:35 -0800
parents
children 1e437c9e5132
comparison
equal deleted inserted replaced
1009:87809b4762c9 1010:f753bada1a46
1 classdef AdvectionRV1D < scheme.Scheme
2 properties
3 grid % Physical grid
4 order % Order accuracy for the approximation
5
6 D % Non-stabalized scheme operator
7 H % Discrete norm
8 Hi % Norm inverse
9 e_l
10 e_r
11 end
12
13 methods
14 function obj = AdvectionRV1D(grid, operator_type, order, waveSpeed)
15 m = grid.size();
16 lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
17 switch operator_type
18 case 'upwind+'
19 ops = sbp.D1Upwind(m, lim, order);
20 D1 = (ops.Dp + ops.Dm)/2;
21 B = ops.e_r*ops.e_r' - ops.e_l*ops.e_l';
22 D2 = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp-ops.HI*(B*spdiag(viscosity)*ops.Dp);
23 DissipationOp = spdiag(abs(waveSpeed))*(ops.Dp-ops.Dm)/2;
24 otherwise
25 error('Other operator types not yet supported', operator_type);
26 end
27 % max(abs()) or just abs()?
28 obj.D = @(viscosity) (-D1 + D2(viscosity) + DissipationOp);
29
30 obj.grid = grid;
31 obj.order = order;
32
33 obj.H = ops.H;
34 obj.Hi = ops.HI;
35 obj.e_l = ops.e_l;
36 obj.e_r = ops.e_r;
37 end
38
39 % Closure functions return the operators applied to the own doamin to close the boundary
40 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
41 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
42 % type is a string specifying the type of boundary condition if there are several.
43 % data is a function returning the data that should be applied at the boundary.
44 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
45 default_arg('type','robin');
46 default_arg('data',0);
47 [e, s] = obj.get_boundary_ops(boundary);
48 switch type
49 case {'D', 'dirichlet'}
50 p = s*obj.Hi*e;
51 closure = p*e';
52 otherwise
53 error('No such boundary condition: type = %s',type);
54 end
55 switch class(data)
56 case 'double'
57 penalty = s*p*data;
58 case 'function_handle'
59 penalty = @(t) s*p*data(t);
60 otherwise
61 error('Wierd data argument!')
62 end
63 end
64
65 % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary.
66 % The right boundary is considered the positive boundary
67 function [e, s] = get_boundary_ops(obj,boundary)
68 switch boundary
69 case 'l'
70 e = obj.e_l;
71 s = -1;
72 case 'r'
73 e = obj.e_r;
74 s = 1;
75 otherwise
76 error('No such boundary: boundary = %s',boundary);
77 end
78 end
79
80 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
81 error('An interface function does not exist yet');
82 end
83
84 function N = size(obj)
85 N = obj.grid.m;
86 end
87 end
88 end