Mercurial > repos > public > sbplib

comparison +scheme/AdvectionRV2D.m @ 1009:87809b4762c9 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Draft implementation of scheme for advection with RV - Implement standard advection scheme. No working support for RV
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 31 Oct 2018 09:31:34 -0700
parents
children
comparison
equal deleted inserted replaced
1008:a6f34de60044 1009:87809b4762c9
1 classdef AdvectionRV2D < scheme.Scheme
2 properties
3 grid % Physical grid
4 order % Order accuracy for the approximation
5
6 D % Non-stabilized scheme operator
7 H % Discrete norm
8 H_inv % Norm inverse
9 halfnorm_inv % Cell array of inverse halfnorm operators
10 e_l % Cell array of left boundary operators
11 e_r % Cell array of right boundary operators
12 d_l % Cell array of left boundary derivative operators
13 d_r % Cell array of right boundary derivative operators
14 waveSpeed
15 end
16
17 methods
18 function obj = AdvectionRV2D(g, operator_type, order, dissipation, waveSpeed)
19 if ~isa(g, 'grid.Cartesian') || g.D() ~= 2
20 error('Grid must be 2d cartesian');
21 end
22
23 obj.grid = g;
24 obj.order = order;
25
26 % Create cell array of 1D operators. For example D1_1d{1} = D1_x, D1_1d{2} = D1_y.
27 [Dp_1d, Dm_1d, H_1d, H_inv_1d, d_l_1d, d_r_1d, e_l_1d, e_r_1d, I, DissipationOp_1d] = ...
28 obj.assemble1DOperators(g, operator_type, order, dissipation);
29
30 %% 2D-operators
31 % D1
32 D1_1d{1} = (Dp_1d{1} + Dp_1d{1})/2;
33 D1_1d{2} = (Dp_1d{2} + Dp_1d{2})/2;
34 D1_2d = obj.extendOperatorTo2D(D1_1d, I);
35 D1 = D1_2d{1} + D1_2d{2};
36 % D2
37
38 Dp_2d = obj.extendOperatorTo2D(Dp_1d, I);
39 Dm_2d = obj.extendOperatorTo2D(Dm_1d, I);
40 D2 = @(viscosity) Dm_2d{1}*spdiag(viscosity)*Dp_2d{1} + Dm_2d{2}*spdiag(viscosity)*Dp_2d{2};
41 % m = g.size();
42 % ind = grid.funcToMatrix(g, 1:g.N());
43 % for i = 1:g.D()
44 % D2_2d{i} = sparse(zeros(g.N()));
45 % end
46 % % x-direction
47 % for i = 1:m(2)
48 % p = ind(:,i);
49 % D2_2d{1}(p,p) = @(viscosity) D2_1d{1}(viscosity(p));
50 % end
51 % % y-direction
52 % for i = 1:m(1)
53 % p = ind(i,:);
54 % D2_2d{2}(p,p) = @(viscosity) D2_1d{2}(viscosity(p));
55 % end
56 % D2 = D2_2d{1} + D2_2d{2};
57
58 obj.d_l = obj.extendOperatorTo2D(d_l_1d, I);
59 obj.d_r = obj.extendOperatorTo2D(d_r_1d, I);
60 obj.e_l = obj.extendOperatorTo2D(e_l_1d, I);
61 obj.e_r = obj.extendOperatorTo2D(e_r_1d, I);
62 obj.H = kron(H_1d{1},H_1d{2});
63 obj.H_inv = kron(H_inv_1d{1},H_inv_1d{2});
64 obj.halfnorm_inv = obj.extendOperatorTo2D(H_inv_1d, I);
65 obj.waveSpeed = waveSpeed;
66
67 % Dissipation operator
68 switch dissipation
69 case 'on'
70 DissOp_2d = obj.extendOperatorTo2D(DissipationOp_1d, I);
71 DissOp = DissOp_2d{1} + DissOp_2d{2};
72 % max(abs()) or just abs()?
73 obj.D = @(v, viscosity) (-waveSpeed.*D1 + D2(viscosity) + abs(waveSpeed).*DissOp)*v;
74 case 'off'
75 obj.D = @(v, viscosity) (-waveSpeed.*D1 + D2(viscosity))*v;
76 end
77 end
78
79 % Closure functions return the operators applied to the own doamin to close the boundary
80 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
81 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
82 % type is a string specifying the type of boundary condition if there are several.
83 % data is a function returning the data that should be applied at the boundary.
84 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
85 default_arg('type','robin');
86 default_arg('data',0);
87 [e, d, halfnorm_inv, i_b, s] = obj.get_boundary_ops(boundary);
88 switch type
89 case {'D', 'dirichlet'}
90 p = s*halfnorm_inv*e;
91 closure = @(v,viscosity) p*(v(i_b));
92 case {'N', 'nuemann'}
93 p = s*halfnorm_inv*e;
94 closure = @(v,viscosity) p*(viscosity(i_b).*d*v);
95 case {'R', 'robin'}
96 p = s*halfnorm_inv*e;
97 closure = @(v, viscosity) p*(obj.waveSpeed(i_b).*v(i_b) - 2*viscosity(i_b).*d*v);
98 otherwise
99 error('No such boundary condition: type = %s',type);
100 end
101 switch class(data)
102 case 'double'
103 penalty = s*p*data;
104 case 'function_handle'
105 penalty = @(t) s*p*data(t);
106 otherwise
107 error('Wierd data argument!')
108 end
109 end
110
111 % Ruturns the boundary ops, half-norm, boundary indices and sign for the boundary specified by the string boundary.
112 % The right boundary for each coordinate direction is considered the positive boundary
113 function [e, d, halfnorm_inv, ind_boundary, s] = get_boundary_ops(obj,boundary)
114 ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N());
115 switch boundary
116 case 'w'
117 e = obj.e_l{1};
118 d = obj.d_l{1};
119 halfnorm_inv = obj.halfnorm_inv{1};
120 ind_boundary = ind(1,:);
121 s = -1;
122 case 'e'
123 e = obj.e_r{1};
124 d = obj.d_r{1};
125 halfnorm_inv = obj.halfnorm_inv{1};
126
127 ind_boundary = ind(end,:);
128 s = 1;
129 case 's'
130 e = obj.e_l{2};
131 d = obj.d_l{2};
132 halfnorm_inv = obj.halfnorm_inv{2};
133 ind_boundary = ind(:,1);
134 s = -1;
135 case 'n'
136 e = obj.e_r{2};
137 d = obj.d_r{2};
138 halfnorm_inv = obj.halfnorm_inv{2};
139 ind_boundary = ind(:,end);
140 s = 1;
141 otherwise
142 error('No such boundary: boundary = %s',boundary);
143 end
144 end
145
146 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
147 error('An interface function does not exist yet');
148 end
149
150 function N = size(obj)
151 N = obj.grid.m;
152 end
153 end
154
155 methods(Static)
156 function [Dp, Dm, H, Hi, d_l, d_r, e_l, e_r, I, DissipationOp] = assemble1DOperators(g, operator_type, order, dissipation)
157 dim = g.D();
158 I = cell(dim,1);
159 D1 = cell(dim,1);
160 D2 = cell(dim,1);
161 H = cell(dim,1);
162 Hi = cell(dim,1);
163 e_l = cell(dim,1);
164 e_r = cell(dim,1);
165 d1_l = cell(dim,1);
166 d1_r = cell(dim,1);
167 DissipationOp = cell(dim,1);
168 for i = 1:dim
169 switch operator_type
170 % case 'narrow'
171 % ops = sbp.D4Variable(g.m(i), g.lim{i}, order);
172 % D1{i} = ops.D1;
173 % D2{i} = ops.D2;
174 % d_l{i} = ops.d1_l';
175 % d_r{i} = ops.d1_r';
176 % if (strcmp(dissipation,'on'))
177 % DissipationOp{i} = -1*sbp.dissipationOperator(g.m(i), order, ops.HI);
178 % end
179 % case 'upwind-'
180 % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order);
181 % D1{i} = (ops.Dp + ops.Dm)/2;
182 % D2{i} = @(viscosity) ops.Dp*spdiag(viscosity)*ops.Dm;
183 % d_l{i} = ops.e_l'*ops.Dm;
184 % d_r{i} = ops.e_r'*ops.Dm;
185 % if (strcmp(dissipation,'on'))
186 % DissipationOp{i} = (ops.Dp-ops.Dm)/2;
187 % end
188 case 'upwind+'
189 ops = sbp.D1Upwind(g.m(i), g.lim{i}, order);
190 Dp{i} = ops.Dp;
191 Dm{i} = ops.Dm;
192 % D1{i} = (ops.Dp + ops.Dm)/2;
193 % D2{i} = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp;
194 d_l{i} = ops.e_l'*ops.Dp;
195 d_r{i} = ops.e_r'*ops.Dp;
196 if (strcmp(dissipation,'on'))
197 DissipationOp{i} = (ops.Dp-ops.Dm)/2;
198 end
199 % case 'upwind+-'
200 % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order);
201 % D1{i} = (ops.Dp + ops.Dm)/2;
202 % D2{i} = @(viscosity) (ops.Dp*spdiag(viscosity)*ops.Dm + ops.Dm*spdiag(viscosity)*ops.Dp)/2;
203 % d_l{i} = ops.e_l'*D1;
204 % d_r{i} = ops.e_r'*D1;
205 % if (strcmp(dissipation,'on'))
206 % DissipationOp{i} = (ops.Dp-ops.Dm)/2;
207 % end
208 otherwise
209 error('Other operator types not yet supported', operator_type);
210 end
211 H{i} = ops.H;
212 Hi{i} = ops.HI;
213 e_l{i} = ops.e_l;
214 e_r{i} = ops.e_r;
215 I{i} = speye(g.m(i));
216 end
217 end
218 function op_2d = extendOperatorTo2D(op, I)
219 op_2d{1} = kr(op{1}, I{2});
220 op_2d{2} = kr(I{1}, op{2});
221 end
222 end
223 end