Mercurial > repos > public > sbplib

comparison +scheme/Utux2d.m @ 1072:6468a5f6ec79 feature/grids/LaplaceSquared

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 12 Feb 2019 17:12:42 +0100
parents 84200bbae101
children 433c89bf19e0
comparison
equal deleted inserted replaced
1071:92cb03e64ca4 1072:6468a5f6ec79
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 H_w, H_e, H_s, H_n % Boundary quadratures
16
17 % Derivatives
18 Dx, Dy
19
20 % Boundary operators
21 e_w, e_e, e_s, e_n
22
23 D % Total discrete operator
24 end
25
26
27 methods
28 function obj = Utux2d(g ,order, opSet, a)
29
30 default_arg('a',1/sqrt(2)*[1, 1]);
31 default_arg('opSet',@sbp.D2Standard);
32
33 assertType(g, 'grid.Cartesian');
34 if iscell(a)
35 a1 = grid.evalOn(g, a{1});
36 a2 = grid.evalOn(g, a{2});
37 a = {spdiag(a1), spdiag(a2)};
38 else
39 a = {a(1), a(2)};
40 end
41
42 m = g.size();
43 m_x = m(1);
44 m_y = m(2);
45 m_tot = g.N();
46
47 xlim = {g.x{1}(1), g.x{1}(end)};
48 ylim = {g.x{2}(1), g.x{2}(end)};
49 obj.grid = g;
50
51 % Operator sets
52 ops_x = opSet(m_x, xlim, order);
53 ops_y = opSet(m_y, ylim, order);
54 Ix = speye(m_x);
55 Iy = speye(m_y);
56
57 % Norms
58 Hx = ops_x.H;
59 Hy = ops_y.H;
60 Hxi = ops_x.HI;
61 Hyi = ops_y.HI;
62
63 obj.H_w = Hy;
64 obj.H_e = Hy;
65 obj.H_s = Hx;
66 obj.H_n = Hx;
67 obj.H_x = Hx;
68 obj.H_y = Hy;
69 obj.H = kron(Hx,Hy);
70 obj.Hi = kron(Hxi,Hyi);
71 obj.Hx = kron(Hx,Iy);
72 obj.Hy = kron(Ix,Hy);
73 obj.Hxi = kron(Hxi,Iy);
74 obj.Hyi = kron(Ix,Hyi);
75
76 % Derivatives
77 Dx = ops_x.D1;
78 Dy = ops_y.D1;
79 obj.Dx = kron(Dx,Iy);
80 obj.Dy = kron(Ix,Dy);
81
82 % Boundary operators
83 obj.e_w = kr(ops_x.e_l, Iy);
84 obj.e_e = kr(ops_x.e_r, Iy);
85 obj.e_s = kr(Ix, ops_y.e_l);
86 obj.e_n = kr(Ix, ops_y.e_r);
87
88 obj.m = m;
89 obj.h = [ops_x.h ops_y.h];
90 obj.order = order;
91 obj.a = a;
92 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
93
94 end
95 % Closure functions return the opertors applied to the own domain to close the boundary
96 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
97 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
98 % type is a string specifying the type of boundary condition if there are several.
99 % data is a function returning the data that should be applied at the boundary.
100 % neighbour_scheme is an instance of Scheme that should be interfaced to.
101 % neighbour_boundary is a string specifying which boundary to interface to.
102 function [closure, penalty] = boundary_condition(obj,boundary,type)
103 default_arg('type','dirichlet');
104
105 sigma = -1; % Scalar penalty parameter
106 switch boundary
107 case {'w','W','west','West'}
108 tau = sigma*obj.a{1}*obj.e_w*obj.H_y;
109 closure = obj.Hi*tau*obj.e_w';
110
111 case {'s','S','south','South'}
112 tau = sigma*obj.a{2}*obj.e_s*obj.H_x;
113 closure = obj.Hi*tau*obj.e_s';
114 end
115 penalty = -obj.Hi*tau;
116
117 end
118
119 % type Struct that specifies the interface coupling.
120 % Fields:
121 % -- couplingType String, type of interface coupling
122 % % Default: 'upwind'. Other: 'centered'
123 % -- interpolation: type of interpolation, default 'none'
124 % -- interpolationDamping: damping on upstream and downstream sides, when using interpolation.
125 % Default {0,0} gives zero damping.
126 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
127
128 defaultType.couplingType = 'upwind';
129 defaultType.interpolation = 'none';
130 defaultType.interpolationDamping = {0,0};
131 default_struct('type', defaultType);
132
133 switch type.interpolation
134 case {'none', ''}
135 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
136 case {'op','OP'}
137 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
138 otherwise
139 error('Unknown type of interpolation: %s ', type.interpolation);
140 end
141 end
142
143 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
144 couplingType = type.couplingType;
145
146 % Get neighbour boundary operator
147 e_neighbour = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
148
149 switch couplingType
150
151 % Upwind coupling (energy dissipation)
152 case 'upwind'
153 sigma_ds = -1; %"Downstream" penalty
154 sigma_us = 0; %"Upstream" penalty
155
156 % Energy-preserving coupling (no energy dissipation)
157 case 'centered'
158 sigma_ds = -1/2; %"Downstream" penalty
159 sigma_us = 1/2; %"Upstream" penalty
160
161 otherwise
162 error(['Interface coupling type ' couplingType ' is not available.'])
163 end
164
165 switch boundary
166 case {'w','W','west','West'}
167 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
168 closure = obj.Hi*tau*obj.e_w';
169 penalty = -obj.Hi*tau*e_neighbour';
170 case {'e','E','east','East'}
171 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
172 closure = obj.Hi*tau*obj.e_e';
173 penalty = -obj.Hi*tau*e_neighbour';
174 case {'s','S','south','South'}
175 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
176 closure = obj.Hi*tau*obj.e_s';
177 penalty = -obj.Hi*tau*e_neighbour';
178 case {'n','N','north','North'}
179 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
180 closure = obj.Hi*tau*obj.e_n';
181 penalty = -obj.Hi*tau*e_neighbour';
182 end
183
184 end
185
186 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
187
188 % User can request special interpolation operators by specifying type.interpOpSet
189 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
190
191 interpOpSet = type.interpOpSet;
192 couplingType = type.couplingType;
193 interpolationDamping = type.interpolationDamping;
194
195 % Get neighbour boundary operator
196 e_neighbour = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
197
198 switch couplingType
199
200 % Upwind coupling (energy dissipation)
201 case 'upwind'
202 sigma_ds = -1; %"Downstream" penalty
203 sigma_us = 0; %"Upstream" penalty
204
205 % Energy-preserving coupling (no energy dissipation)
206 case 'centered'
207 sigma_ds = -1/2; %"Downstream" penalty
208 sigma_us = 1/2; %"Upstream" penalty
209
210 otherwise
211 error(['Interface coupling type ' couplingType ' is not available.'])
212 end
213
214 int_damp_us = interpolationDamping{1};
215 int_damp_ds = interpolationDamping{2};
216
217 % u denotes the solution in the own domain
218 % v denotes the solution in the neighbour domain
219 % Find the number of grid points along the interface
220 switch boundary
221 case {'w','e'}
222 m_u = obj.m(2);
223 case {'s','n'}
224 m_u = obj.m(1);
225 end
226 m_v = size(e_neighbour, 2);
227
228 % Build interpolation operators
229 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
230 Iu2v = intOps.Iu2v;
231 Iv2u = intOps.Iv2u;
232
233 I_local2neighbour_ds = intOps.Iu2v.bad;
234 I_local2neighbour_us = intOps.Iu2v.good;
235 I_neighbour2local_ds = intOps.Iv2u.good;
236 I_neighbour2local_us = intOps.Iv2u.bad;
237
238 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
239 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
240
241
242 switch boundary
243 case {'w','W','west','West'}
244 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
245 closure = obj.Hi*tau*obj.e_w';
246 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
247
248 beta = int_damp_ds*obj.a{1}...
249 *obj.e_w*obj.H_y;
250 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_w' - obj.Hi*beta*obj.e_w';
251 case {'e','E','east','East'}
252 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
253 closure = obj.Hi*tau*obj.e_e';
254 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
255
256 beta = int_damp_us*obj.a{1}...
257 *obj.e_e*obj.H_y;
258 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_e' - obj.Hi*beta*obj.e_e';
259 case {'s','S','south','South'}
260 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
261 closure = obj.Hi*tau*obj.e_s';
262 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
263
264 beta = int_damp_ds*obj.a{2}...
265 *obj.e_s*obj.H_x;
266 closure = closure + obj.Hi*beta*I_back_forth_ds*obj.e_s' - obj.Hi*beta*obj.e_s';
267 case {'n','N','north','North'}
268 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
269 closure = obj.Hi*tau*obj.e_n';
270 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
271
272 beta = int_damp_us*obj.a{2}...
273 *obj.e_n*obj.H_x;
274 closure = closure + obj.Hi*beta*I_back_forth_us*obj.e_n' - obj.Hi*beta*obj.e_n';
275 end
276
277
278 end
279
280 % Returns the boundary operator op for the boundary specified by the string boundary.
281 % op -- string
282 % boundary -- string
283 function o = getBoundaryOperator(obj, op, boundary)
284 assertIsMember(op, {'e'})
285 assertIsMember(boundary, {'w', 'e', 's', 'n'})
286
287 o = obj.([op, '_', boundary]);
288 end
289
290 % Returns square boundary quadrature matrix, of dimension
291 % corresponding to the number of boundary points
292 %
293 % boundary -- string
294 function H_b = getBoundaryQuadrature(obj, boundary)
295 assertIsMember(boundary, {'w', 'e', 's', 'n'})
296
297 H_b = obj.(['H_', boundary]);
298 end
299
300 function N = size(obj)
301 N = obj.m;
302 end
303
304 end
305 end