Mercurial > repos > public > sbplib

comparison +scheme/Utux2d.m @ 950:cab047de7f5d feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename *2D schemes to *2d
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 06 Dec 2018 10:32:02 +0100
parents +scheme/Utux2D.m@3dd7f87c9a1b
children 78db023a7fe3 f029b97dbc72 c12b84fe9b00
comparison
equal deleted inserted replaced
949:6d2167719557 950:cab047de7f5d
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 v0 % Initial data
8
9 a % Wave speed a = [a1, a2];
10 % Can either be a constant vector or a cell array of function handles.
11
12 H % Discrete norm
13 H_x, H_y % Norms in the x and y directions
14 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
15
16 % Derivatives
17 Dx, Dy
18
19 % Boundary operators
20 e_w, e_e, e_s, e_n
21
22 D % Total discrete operator
23 end
24
25
26 methods
27 function obj = Utux2d(g ,order, opSet, a)
28
29 default_arg('a',1/sqrt(2)*[1, 1]);
30 default_arg('opSet',@sbp.D2Standard);
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 Dx = ops_x.D1;
73 Dy = ops_y.D1;
74 obj.Dx = kron(Dx,Iy);
75 obj.Dy = kron(Ix,Dy);
76
77 % Boundary operators
78 obj.e_w = kr(ops_x.e_l, Iy);
79 obj.e_e = kr(ops_x.e_r, Iy);
80 obj.e_s = kr(Ix, ops_y.e_l);
81 obj.e_n = kr(Ix, ops_y.e_r);
82
83 obj.m = m;
84 obj.h = [ops_x.h ops_y.h];
85 obj.order = order;
86 obj.a = a;
87 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
88
89 end
90 % Closure functions return the opertors applied to the own domain to close the boundary
91 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
92 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
93 % type is a string specifying the type of boundary condition if there are several.
94 % data is a function returning the data that should be applied at the boundary.
95 % neighbour_scheme is an instance of Scheme that should be interfaced to.
96 % neighbour_boundary is a string specifying which boundary to interface to.
97 function [closure, penalty] = boundary_condition(obj,boundary,type)
98 default_arg('type','dirichlet');
99
100 sigma = -1; % Scalar penalty parameter
101 switch boundary
102 case {'w','W','west','West'}
103 tau = sigma*obj.a{1}*obj.e_w*obj.H_y;
104 closure = obj.Hi*tau*obj.e_w';
105
106 case {'s','S','south','South'}
107 tau = sigma*obj.a{2}*obj.e_s*obj.H_x;
108 closure = obj.Hi*tau*obj.e_s';
109 end
110 penalty = -obj.Hi*tau;
111
112 end
113
114 % type Struct that specifies the interface coupling.
115 % Fields:
116 % -- couplingType String, type of interface coupling
117 % % Default: 'upwind'. Other: 'centered'
118 % -- interpolation: type of interpolation, default 'none'
119 % -- interpolationDamping: damping on upstream and downstream sides, when using interpolation.
120 % Default {0,0} gives zero damping.
121 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
122
123 defaultType.couplingType = 'upwind';
124 defaultType.interpolation = 'none';
125 defaultType.interpolationDamping = {0,0};
126 default_struct('type', defaultType);
127
128 switch type.interpolation
129 case {'none', ''}
130 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
131 case {'op','OP'}
132 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
133 otherwise
134 error('Unknown type of interpolation: %s ', type.interpolation);
135 end
136 end
137
138 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
139 couplingType = type.couplingType;
140
141 % Get neighbour boundary operator
142 switch neighbour_boundary
143 case {'e','E','east','East'}
144 e_neighbour = neighbour_scheme.e_e;
145 case {'w','W','west','West'}
146 e_neighbour = neighbour_scheme.e_w;
147 case {'n','N','north','North'}
148 e_neighbour = neighbour_scheme.e_n;
149 case {'s','S','south','South'}
150 e_neighbour = neighbour_scheme.e_s;
151 end
152
153 switch couplingType
154
155 % Upwind coupling (energy dissipation)
156 case 'upwind'
157 sigma_ds = -1; %"Downstream" penalty
158 sigma_us = 0; %"Upstream" penalty
159
160 % Energy-preserving coupling (no energy dissipation)
161 case 'centered'
162 sigma_ds = -1/2; %"Downstream" penalty
163 sigma_us = 1/2; %"Upstream" penalty
164
165 otherwise
166 error(['Interface coupling type ' couplingType ' is not available.'])
167 end
168
169 switch boundary
170 case {'w','W','west','West'}
171 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
172 closure = obj.Hi*tau*obj.e_w';
173 penalty = -obj.Hi*tau*e_neighbour';
174 case {'e','E','east','East'}
175 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
176 closure = obj.Hi*tau*obj.e_e';
177 penalty = -obj.Hi*tau*e_neighbour';
178 case {'s','S','south','South'}
179 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
180 closure = obj.Hi*tau*obj.e_s';
181 penalty = -obj.Hi*tau*e_neighbour';
182 case {'n','N','north','North'}
183 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
184 closure = obj.Hi*tau*obj.e_n';
185 penalty = -obj.Hi*tau*e_neighbour';
186 end
187
188 end
189
190 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
191
192 % User can request special interpolation operators by specifying type.interpOpSet
193 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
194
195 interpOpSet = type.interpOpSet;
196 couplingType = type.couplingType;
197 interpolationDamping = type.interpolationDamping;
198
199 % Get neighbour boundary operator
200 switch neighbour_boundary
201 case {'e','E','east','East'}
202 e_neighbour = neighbour_scheme.e_e;
203 case {'w','W','west','West'}
204 e_neighbour = neighbour_scheme.e_w;
205 case {'n','N','north','North'}
206 e_neighbour = neighbour_scheme.e_n;
207 case {'s','S','south','South'}
208 e_neighbour = neighbour_scheme.e_s;
209 end
210
211 switch couplingType
212
213 % Upwind coupling (energy dissipation)
214 case 'upwind'
215 sigma_ds = -1; %"Downstream" penalty
216 sigma_us = 0; %"Upstream" penalty
217
218 % Energy-preserving coupling (no energy dissipation)
219 case 'centered'
220 sigma_ds = -1/2; %"Downstream" penalty
221 sigma_us = 1/2; %"Upstream" penalty
222
223 otherwise
224 error(['Interface coupling type ' couplingType ' is not available.'])
225 end
226
227 int_damp_us = interpolationDamping{1};
228 int_damp_ds = interpolationDamping{2};
229
230 % u denotes the solution in the own domain
231 % v denotes the solution in the neighbour domain
232 % Find the number of grid points along the interface
233 switch boundary
234 case {'w','e'}
235 m_u = obj.m(2);
236 case {'s','n'}
237 m_u = obj.m(1);
238 end
239 m_v = size(e_neighbour, 2);
240
241 % Build interpolation operators
242 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
243 Iu2v = intOps.Iu2v;
244 Iv2u = intOps.Iv2u;
245
246 I_local2neighbour_ds = intOps.Iu2v.bad;
247 I_local2neighbour_us = intOps.Iu2v.good;
248 I_neighbour2local_ds = intOps.Iv2u.good;
249 I_neighbour2local_us = intOps.Iv2u.bad;
250
251 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
252 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
253
254
255 switch boundary
256 case {'w','W','west','West'}
257 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
258 closure = obj.Hi*tau*obj.e_w';
259 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
260
261 beta = int_damp_ds*obj.a{1}...
262 *obj.e_w*obj.H_y;
263 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_w' - obj.Hi*beta*obj.e_w';
264 case {'e','E','east','East'}
265 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
266 closure = obj.Hi*tau*obj.e_e';
267 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
268
269 beta = int_damp_us*obj.a{1}...
270 *obj.e_e*obj.H_y;
271 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_e' - obj.Hi*beta*obj.e_e';
272 case {'s','S','south','South'}
273 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
274 closure = obj.Hi*tau*obj.e_s';
275 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
276
277 beta = int_damp_ds*obj.a{2}...
278 *obj.e_s*obj.H_x;
279 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_s' - obj.Hi*beta*obj.e_s';
280 case {'n','N','north','North'}
281 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
282 closure = obj.Hi*tau*obj.e_n';
283 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
284
285 beta = int_damp_us*obj.a{2}...
286 *obj.e_n*obj.H_x;
287 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_n' - obj.Hi*beta*obj.e_n';
288 end
289
290
291 end
292
293 function N = size(obj)
294 N = obj.m;
295 end
296
297 end
298
299 methods(Static)
300 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
301 % and bound_v of scheme schm_v.
302 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
303 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
304 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
305 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
306 end
307 end
308 end