Mercurial > repos > public > sbplib

comparison +scheme/Utux2d.m @ 1033:037f203b9bf5 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with branch feature/advectioRV to utilize the +rv package
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:44:12 +0100
parents ac80bedc8df7
children 8a9393084b30
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1033:037f203b9bf5
1 classdef Utux2d < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 grid % Grid
6 order % Order accuracy for the approximation
7
8 a % Wave speed a = [a1, a2];
9 % Can either be a constant vector or a cell array of function handles.
10
11 H % Discrete norm
12 H_x, H_y % Norms in the x and y directions
13 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
14
15 % Derivatives
16 Dx, Dy
17
18 % Boundary operators
19 e_w, e_e, e_s, e_n
20
21 D % Total discrete operator
22 end
23
24
25 methods
26 function obj = Utux2d(g ,order, opSet, a, fluxSplitting)
27
28 default_arg('a',1/sqrt(2)*[1, 1]);
29 default_arg('opSet',@sbp.D2Standard);
30 default_arg('fluxSplitting',[]);
31
32 assertType(g, 'grid.Cartesian');
33 if iscell(a)
34 a1 = grid.evalOn(g, a{1});
35 a2 = grid.evalOn(g, a{2});
36 a = {spdiag(a1), spdiag(a2)};
37 else
38 a = {a(1), a(2)};
39 end
40
41 m = g.size();
42 m_x = m(1);
43 m_y = m(2);
44 m_tot = g.N();
45
46 xlim = {g.x{1}(1), g.x{1}(end)};
47 ylim = {g.x{2}(1), g.x{2}(end)};
48 obj.grid = g;
49
50 % Operator sets
51 ops_x = opSet(m_x, xlim, order);
52 ops_y = opSet(m_y, ylim, order);
53 Ix = speye(m_x);
54 Iy = speye(m_y);
55
56 % Norms
57 Hx = ops_x.H;
58 Hy = ops_y.H;
59 Hxi = ops_x.HI;
60 Hyi = ops_y.HI;
61
62 obj.H_x = Hx;
63 obj.H_y = Hy;
64 obj.H = kron(Hx,Hy);
65 obj.Hi = kron(Hxi,Hyi);
66 obj.Hx = kron(Hx,Iy);
67 obj.Hy = kron(Ix,Hy);
68 obj.Hxi = kron(Hxi,Iy);
69 obj.Hyi = kron(Ix,Hyi);
70
71 % Derivatives
72 if (isequal(opSet,@sbp.D1Upwind))
73 Dx = (ops_x.Dp + ops_x.Dm)/2;
74 Dy = (ops_y.Dp + ops_y.Dm)/2;
75 obj.Dx = kron(Dx,Iy);
76 obj.Dy = kron(Ix,Dy);
77 DissOpx = (ops_x.Dm - ops_x.Dp)/2;
78 DissOpy = (ops_y.Dm - ops_y.Dp)/2;
79 DissOpx = kron(DissOpx,Iy);
80 DissOpy = kron(Ix,DissOpy);
81
82 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy + fluxSplitting{1}*DissOpx + fluxSplitting{2}*DissOpy);
83 else
84 Dx = ops_x.D1;
85 Dy = ops_y.D1;
86 obj.Dx = kron(Dx,Iy);
87 obj.Dy = kron(Ix,Dy);
88
89 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
90 end
91
92 % Boundary operators
93 obj.e_w = kr(ops_x.e_l, Iy);
94 obj.e_e = kr(ops_x.e_r, Iy);
95 obj.e_s = kr(Ix, ops_y.e_l);
96 obj.e_n = kr(Ix, ops_y.e_r);
97
98 obj.m = m;
99 obj.h = [ops_x.h ops_y.h];
100 obj.order = order;
101 obj.a = a;
102 end
103 % Closure functions return the opertors applied to the own domain to close the boundary
104 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
105 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
106 % type is a string specifying the type of boundary condition if there are several. %TBD Remove type here? Only dirichlet applicable?
107 % data is a function returning the data that should be applied at the boundary.
108 % neighbour_scheme is an instance of Scheme that should be interfaced to.
109 % neighbour_boundary is a string specifying which boundary to interface to.
110 function [closure, penalty] = boundary_condition(obj,boundary,type)
111 default_arg('type','dirichlet');
112 sigma_left = -1; % Scalar penalty parameter for left boundaries (West/South)
113 sigma_right = 1; % Scalar penalty parameter for right boundaries (East/North)
114 switch boundary
115 % Can only specify boundary condition where there is inflow
116 % Extract the postivie resp. negative part of a, for the left
117 % resp. right boundaries, and set other values of a to zero.
118 % Then the closure will effectively only contribute to inflow boundaries
119 case {'w','W','west','West'}
120 a_inflow = obj.a{1};
121 a_inflow(a_inflow < 0) = 0;
122 tau = sigma_left*a_inflow*obj.e_w*obj.H_y;
123 closure = obj.Hi*tau*obj.e_w';
124 case {'e','E','east','East'}
125 a_inflow = obj.a{1};
126 a_inflow(a_inflow > 0) = 0;
127 tau = sigma_right*a_inflow*obj.e_e*obj.H_y;
128 closure = obj.Hi*tau*obj.e_e';
129 case {'s','S','south','South'}
130 a_inflow = obj.a{2};
131 a_inflow(a_inflow < 0) = 0;
132 tau = sigma_left*a_inflow*obj.e_s*obj.H_x;
133 closure = obj.Hi*tau*obj.e_s';
134 case {'n','N','north','North'}
135 a_inflow = obj.a{2};
136 a_inflow(a_inflow > 0) = 0;
137 tau = sigma_right*a_inflow*obj.e_n*obj.H_x;
138 closure = obj.Hi*tau*obj.e_n';
139 end
140 penalty = -obj.Hi*tau;
141 end
142
143 % type Struct that specifies the interface coupling.
144 % Fields:
145 % -- couplingType String, type of interface coupling
146 % % Default: 'upwind'. Other: 'centered'
147 % -- interpolation: type of interpolation, default 'none'
148 % -- interpolationDamping: damping on upstream and downstream sides, when using interpolation.
149 % Default {0,0} gives zero damping.
150 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
151
152 defaultType.couplingType = 'upwind';
153 defaultType.interpolation = 'none';
154 defaultType.interpolationDamping = {0,0};
155 default_struct('type', defaultType);
156
157 switch type.interpolation
158 case {'none', ''}
159 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
160 case {'op','OP'}
161 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
162 otherwise
163 error('Unknown type of interpolation: %s ', type.interpolation);
164 end
165 end
166
167 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
168 couplingType = type.couplingType;
169
170 % Get neighbour boundary operator
171 switch neighbour_boundary
172 case {'e','E','east','East'}
173 e_neighbour = neighbour_scheme.e_e;
174 case {'w','W','west','West'}
175 e_neighbour = neighbour_scheme.e_w;
176 case {'n','N','north','North'}
177 e_neighbour = neighbour_scheme.e_n;
178 case {'s','S','south','South'}
179 e_neighbour = neighbour_scheme.e_s;
180 end
181
182 switch couplingType
183
184 % Upwind coupling (energy dissipation)
185 case 'upwind'
186 sigma_ds = -1; %"Downstream" penalty
187 sigma_us = 0; %"Upstream" penalty
188
189 % Energy-preserving coupling (no energy dissipation)
190 case 'centered'
191 sigma_ds = -1/2; %"Downstream" penalty
192 sigma_us = 1/2; %"Upstream" penalty
193
194 otherwise
195 error(['Interface coupling type ' couplingType ' is not available.'])
196 end
197
198 switch boundary
199 case {'w','W','west','West'}
200 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
201 closure = obj.Hi*tau*obj.e_w';
202 penalty = -obj.Hi*tau*e_neighbour';
203 case {'e','E','east','East'}
204 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
205 closure = obj.Hi*tau*obj.e_e';
206 penalty = -obj.Hi*tau*e_neighbour';
207 case {'s','S','south','South'}
208 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
209 closure = obj.Hi*tau*obj.e_s';
210 penalty = -obj.Hi*tau*e_neighbour';
211 case {'n','N','north','North'}
212 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
213 closure = obj.Hi*tau*obj.e_n';
214 penalty = -obj.Hi*tau*e_neighbour';
215 end
216
217 end
218
219 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
220
221 % User can request special interpolation operators by specifying type.interpOpSet
222 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
223
224 interpOpSet = type.interpOpSet;
225 couplingType = type.couplingType;
226 interpolationDamping = type.interpolationDamping;
227
228 % Get neighbour boundary operator
229 switch neighbour_boundary
230 case {'e','E','east','East'}
231 e_neighbour = neighbour_scheme.e_e;
232 case {'w','W','west','West'}
233 e_neighbour = neighbour_scheme.e_w;
234 case {'n','N','north','North'}
235 e_neighbour = neighbour_scheme.e_n;
236 case {'s','S','south','South'}
237 e_neighbour = neighbour_scheme.e_s;
238 end
239
240 switch couplingType
241
242 % Upwind coupling (energy dissipation)
243 case 'upwind'
244 sigma_ds = -1; %"Downstream" penalty
245 sigma_us = 0; %"Upstream" penalty
246
247 % Energy-preserving coupling (no energy dissipation)
248 case 'centered'
249 sigma_ds = -1/2; %"Downstream" penalty
250 sigma_us = 1/2; %"Upstream" penalty
251
252 otherwise
253 error(['Interface coupling type ' couplingType ' is not available.'])
254 end
255
256 int_damp_us = interpolationDamping{1};
257 int_damp_ds = interpolationDamping{2};
258
259 % u denotes the solution in the own domain
260 % v denotes the solution in the neighbour domain
261 % Find the number of grid points along the interface
262 switch boundary
263 case {'w','e'}
264 m_u = obj.m(2);
265 case {'s','n'}
266 m_u = obj.m(1);
267 end
268 m_v = size(e_neighbour, 2);
269
270 % Build interpolation operators
271 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
272 Iu2v = intOps.Iu2v;
273 Iv2u = intOps.Iv2u;
274
275 I_local2neighbour_ds = intOps.Iu2v.bad;
276 I_local2neighbour_us = intOps.Iu2v.good;
277 I_neighbour2local_ds = intOps.Iv2u.good;
278 I_neighbour2local_us = intOps.Iv2u.bad;
279
280 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
281 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
282
283
284 switch boundary
285 case {'w','W','west','West'}
286 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
287 closure = obj.Hi*tau*obj.e_w';
288 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
289
290 beta = int_damp_ds*obj.a{1}...
291 *obj.e_w*obj.H_y;
292 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_w' - obj.Hi*beta*obj.e_w';
293 case {'e','E','east','East'}
294 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
295 closure = obj.Hi*tau*obj.e_e';
296 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
297
298 beta = int_damp_us*obj.a{1}...
299 *obj.e_e*obj.H_y;
300 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_e' - obj.Hi*beta*obj.e_e';
301 case {'s','S','south','South'}
302 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
303 closure = obj.Hi*tau*obj.e_s';
304 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
305
306 beta = int_damp_ds*obj.a{2}...
307 *obj.e_s*obj.H_x;
308 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_s' - obj.Hi*beta*obj.e_s';
309 case {'n','N','north','North'}
310 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
311 closure = obj.Hi*tau*obj.e_n';
312 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
313
314 beta = int_damp_us*obj.a{2}...
315 *obj.e_n*obj.H_x;
316 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_n' - obj.Hi*beta*obj.e_n';
317 end
318
319
320 end
321
322 function N = size(obj)
323 N = obj.m;
324 end
325
326 end
327
328 methods(Static)
329 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
330 % and bound_v of scheme schm_v.
331 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
332 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
333 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
334 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
335 end
336 end
337 end