Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearVirtaMin.m @ 1136:eee71789f13b feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add LaplaceCurvilinear schemes where the minimum check will be implemented. The Virta scheme will be used for comparison only.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 10 Jun 2019 10:43:12 +0200
parents
children 2ff1f366e64a
comparison
equal deleted inserted replaced
1135:a0c6f060a105 1136:eee71789f13b
1 classdef LaplaceCurvilinearVirtaMin < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5
6 grid
7
8 order % Order accuracy for the approximation
9
10 a,b % Parameters of the operator
11
12
13 % Inner products and operators for physical coordinates
14 D % Laplace operator
15 H, Hi % Inner product
16 e_w, e_e, e_s, e_n
17 d_w, d_e, d_s, d_n % Normal derivatives at the boundary
18 H_w, H_e, H_s, H_n % Boundary inner products
19 Dx, Dy % Physical derivatives
20 M % Gradient inner product
21
22 % Metric coefficients
23 J, Ji
24 a11, a12, a22
25 x_u
26 x_v
27 y_u
28 y_v
29
30 % Inner product and operators for logical coordinates
31 H_u, H_v % Norms in the x and y directions
32 Hi_u, Hi_v
33 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
34 Hiu, Hiv
35 du_w, dv_w
36 du_e, dv_e
37 du_s, dv_s
38 du_n, dv_n
39 gamm_u, gamm_v
40 lambda
41
42 end
43
44 methods
45 % Implements a*div(b*grad(u)) as a SBP scheme
46 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
47
48 function obj = LaplaceCurvilinearVirtaMin(g ,order, a, b, opSet)
49 default_arg('opSet',@sbp.D2Variable);
50 default_arg('a', 1);
51 default_arg('b', 1);
52
53 if b ~=1
54 error('Not implemented yet')
55 end
56
57 % assert(isa(g, 'grid.Curvilinear'))
58 if isa(a, 'function_handle')
59 a = grid.evalOn(g, a);
60 a = spdiag(a);
61 end
62
63 m = g.size();
64 m_u = m(1);
65 m_v = m(2);
66 m_tot = g.N();
67
68 h = g.scaling();
69 h_u = h(1);
70 h_v = h(2);
71
72
73 % 1D operators
74 ops_u = opSet(m_u, {0, 1}, order);
75 ops_v = opSet(m_v, {0, 1}, order);
76
77 I_u = speye(m_u);
78 I_v = speye(m_v);
79
80 D1_u = ops_u.D1;
81 D2_u = ops_u.D2;
82 H_u = ops_u.H;
83 Hi_u = ops_u.HI;
84 e_l_u = ops_u.e_l;
85 e_r_u = ops_u.e_r;
86 d1_l_u = ops_u.d1_l;
87 d1_r_u = ops_u.d1_r;
88
89 D1_v = ops_v.D1;
90 D2_v = ops_v.D2;
91 H_v = ops_v.H;
92 Hi_v = ops_v.HI;
93 e_l_v = ops_v.e_l;
94 e_r_v = ops_v.e_r;
95 d1_l_v = ops_v.d1_l;
96 d1_r_v = ops_v.d1_r;
97
98
99 % Logical operators
100 Du = kr(D1_u,I_v);
101 Dv = kr(I_u,D1_v);
102 obj.Hu = kr(H_u,I_v);
103 obj.Hv = kr(I_u,H_v);
104 obj.Hiu = kr(Hi_u,I_v);
105 obj.Hiv = kr(I_u,Hi_v);
106
107 e_w = kr(e_l_u,I_v);
108 e_e = kr(e_r_u,I_v);
109 e_s = kr(I_u,e_l_v);
110 e_n = kr(I_u,e_r_v);
111 obj.du_w = kr(d1_l_u,I_v);
112 obj.dv_w = (e_w'*Dv)';
113 obj.du_e = kr(d1_r_u,I_v);
114 obj.dv_e = (e_e'*Dv)';
115 obj.du_s = (e_s'*Du)';
116 obj.dv_s = kr(I_u,d1_l_v);
117 obj.du_n = (e_n'*Du)';
118 obj.dv_n = kr(I_u,d1_r_v);
119
120
121 % Metric coefficients
122 coords = g.points();
123 x = coords(:,1);
124 y = coords(:,2);
125
126 x_u = Du*x;
127 x_v = Dv*x;
128 y_u = Du*y;
129 y_v = Dv*y;
130
131 J = x_u.*y_v - x_v.*y_u;
132 a11 = 1./J .* (x_v.^2 + y_v.^2);
133 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
134 a22 = 1./J .* (x_u.^2 + y_u.^2);
135 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
136
137 obj.x_u = x_u;
138 obj.x_v = x_v;
139 obj.y_u = y_u;
140 obj.y_v = y_v;
141
142
143 % Assemble full operators
144 L_12 = spdiag(a12);
145 Duv = Du*L_12*Dv;
146 Dvu = Dv*L_12*Du;
147
148 Duu = sparse(m_tot);
149 Dvv = sparse(m_tot);
150 ind = grid.funcToMatrix(g, 1:m_tot);
151
152 for i = 1:m_v
153 D = D2_u(a11(ind(:,i)));
154 p = ind(:,i);
155 Duu(p,p) = D;
156 end
157
158 for i = 1:m_u
159 D = D2_v(a22(ind(i,:)));
160 p = ind(i,:);
161 Dvv(p,p) = D;
162 end
163
164
165 % Physical operators
166 obj.J = spdiag(J);
167 obj.Ji = spdiag(1./J);
168
169 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv);
170 obj.H = obj.J*kr(H_u,H_v);
171 obj.Hi = obj.Ji*kr(Hi_u,Hi_v);
172
173 obj.e_w = e_w;
174 obj.e_e = e_e;
175 obj.e_s = e_s;
176 obj.e_n = e_n;
177
178 %% normal derivatives
179 I_w = ind(1,:);
180 I_e = ind(end,:);
181 I_s = ind(:,1);
182 I_n = ind(:,end);
183
184 a11_w = spdiag(a11(I_w));
185 a12_w = spdiag(a12(I_w));
186 a11_e = spdiag(a11(I_e));
187 a12_e = spdiag(a12(I_e));
188 a22_s = spdiag(a22(I_s));
189 a12_s = spdiag(a12(I_s));
190 a22_n = spdiag(a22(I_n));
191 a12_n = spdiag(a12(I_n));
192
193 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
194 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
195 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
196 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
197
198 obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))';
199 obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))';
200 obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))';
201 obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))';
202
203 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
204 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
205
206 %% Boundary inner products
207 obj.H_w = H_v*spdiag(s_w);
208 obj.H_e = H_v*spdiag(s_e);
209 obj.H_s = H_u*spdiag(s_s);
210 obj.H_n = H_u*spdiag(s_n);
211
212 % Misc.
213 obj.m = m;
214 obj.h = [h_u h_v];
215 obj.order = order;
216 obj.grid = g;
217
218 obj.a = a;
219 obj.b = b;
220 obj.a11 = a11;
221 obj.a12 = a12;
222 obj.a22 = a22;
223 obj.lambda = lambda;
224
225 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
226 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
227 end
228
229
230 % Closure functions return the opertors applied to the own doamin to close the boundary
231 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
232 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
233 % type is a string specifying the type of boundary condition if there are several.
234 % data is a function returning the data that should be applied at the boundary.
235 % neighbour_scheme is an instance of Scheme that should be interfaced to.
236 % neighbour_boundary is a string specifying which boundary to interface to.
237 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
238 default_arg('type','neumann');
239 default_arg('parameter', []);
240
241 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
242 H_b = obj.getBoundaryQuadrature(boundary);
243 gamm = obj.getBoundaryBorrowing(boundary);
244
245 switch type
246 % Dirichlet boundary condition
247 case {'D','d','dirichlet'}
248 tuning = 1.0;
249
250 b1 = gamm*obj.lambda./obj.a11.^2;
251 b2 = gamm*obj.lambda./obj.a22.^2;
252
253 tau1 = tuning * spdiag(-1./b1 - 1./b2);
254 tau2 = 1;
255
256 tau = (tau1*e + tau2*d)*H_b;
257
258 closure = obj.a*obj.Hi*tau*e';
259 penalty = -obj.a*obj.Hi*tau;
260
261
262 % Neumann boundary condition
263 case {'N','n','neumann'}
264 tau1 = -1;
265 tau2 = 0;
266 tau = (tau1*e + tau2*d)*H_b;
267
268 closure = obj.a*obj.Hi*tau*d';
269 penalty = -obj.a*obj.Hi*tau;
270
271
272 % Unknown, boundary condition
273 otherwise
274 error('No such boundary condition: type = %s',type);
275 end
276 end
277
278 % type Struct that specifies the interface coupling.
279 % Fields:
280 % -- tuning: penalty strength, defaults to 1.2
281 % -- interpolation: type of interpolation, default 'none'
282 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
283
284 defaultType.tuning = 1.2;
285 defaultType.interpolation = 'none';
286 default_struct('type', defaultType);
287
288 switch type.interpolation
289 case {'none', ''}
290 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
291 case {'op','OP'}
292 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
293 otherwise
294 error('Unknown type of interpolation: %s ', type.interpolation);
295 end
296 end
297
298 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
299 tuning = type.tuning;
300
301 % u denotes the solution in the own domain
302 % v denotes the solution in the neighbour domain
303 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
304 H_b_u = obj.getBoundaryQuadrature(boundary);
305 I_u = obj.getBoundaryIndices(boundary);
306 gamm_u = obj.getBoundaryBorrowing(boundary);
307
308 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
309 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
310 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
311 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
312
313 u = obj;
314 v = neighbour_scheme;
315
316 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
317 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
318 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
319 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
320
321 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
322 tau1 = tuning * spdiag(tau1);
323 tau2 = 1/2;
324
325 sig1 = -1/2;
326 sig2 = 0;
327
328 tau = (e_u*tau1 + tau2*d_u)*H_b_u;
329 sig = (sig1*e_u + sig2*d_u)*H_b_u;
330
331 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
332 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
333 end
334
335 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
336
337 % TODO: Make this work for curvilinear grids
338 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
339
340 % User can request special interpolation operators by specifying type.interpOpSet
341 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
342 interpOpSet = type.interpOpSet;
343 tuning = type.tuning;
344
345
346 % u denotes the solution in the own domain
347 % v denotes the solution in the neighbour domain
348 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
349 H_b_u = obj.getBoundaryQuadrature(boundary);
350 I_u = obj.getBoundaryIndices(boundary);
351 gamm_u = obj.getBoundaryBorrowing(boundary);
352
353 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
354 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
355 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
356 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
357
358
359 % Find the number of grid points along the interface
360 m_u = size(e_u, 2);
361 m_v = size(e_v, 2);
362
363 Hi = obj.Hi;
364 a = obj.a;
365
366 u = obj;
367 v = neighbour_scheme;
368
369 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
370 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
371 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
372 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
373
374 tau_u = -1./(4*b1_u) -1./(4*b2_u);
375 tau_v = -1./(4*b1_v) -1./(4*b2_v);
376
377 tau_u = tuning * spdiag(tau_u);
378 tau_v = tuning * spdiag(tau_v);
379 beta_u = tau_v;
380
381 % Build interpolation operators
382 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
383 Iu2v = intOps.Iu2v;
384 Iv2u = intOps.Iv2u;
385
386 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
387 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
388 a*1/2*Hi*d_u*H_b_u*e_u' + ...
389 -a*1/2*Hi*e_u*H_b_u*d_u';
390
391 penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ...
392 -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ...
393 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
394 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
395
396 end
397
398 % Returns the boundary operator op for the boundary specified by the string boundary.
399 % op -- string or a cell array of strings
400 % boundary -- string
401 function varargout = getBoundaryOperator(obj, op, boundary)
402
403 if ~iscell(op)
404 op = {op};
405 end
406
407 for i = 1:numel(op)
408 switch op{i}
409 case 'e'
410 switch boundary
411 case 'w'
412 e = obj.e_w;
413 case 'e'
414 e = obj.e_e;
415 case 's'
416 e = obj.e_s;
417 case 'n'
418 e = obj.e_n;
419 otherwise
420 error('No such boundary: boundary = %s',boundary);
421 end
422 varargout{i} = e;
423
424 case 'd'
425 switch boundary
426 case 'w'
427 d = obj.d_w;
428 case 'e'
429 d = obj.d_e;
430 case 's'
431 d = obj.d_s;
432 case 'n'
433 d = obj.d_n;
434 otherwise
435 error('No such boundary: boundary = %s',boundary);
436 end
437 varargout{i} = d;
438 end
439 end
440 end
441
442 % Returns square boundary quadrature matrix, of dimension
443 % corresponding to the number of boundary points
444 %
445 % boundary -- string
446 function H_b = getBoundaryQuadrature(obj, boundary)
447
448 switch boundary
449 case 'w'
450 H_b = obj.H_w;
451 case 'e'
452 H_b = obj.H_e;
453 case 's'
454 H_b = obj.H_s;
455 case 'n'
456 H_b = obj.H_n;
457 otherwise
458 error('No such boundary: boundary = %s',boundary);
459 end
460 end
461
462 % Returns the indices of the boundary points in the grid matrix
463 % boundary -- string
464 function I = getBoundaryIndices(obj, boundary)
465 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
466 switch boundary
467 case 'w'
468 I = ind(1,:);
469 case 'e'
470 I = ind(end,:);
471 case 's'
472 I = ind(:,1)';
473 case 'n'
474 I = ind(:,end)';
475 otherwise
476 error('No such boundary: boundary = %s',boundary);
477 end
478 end
479
480 % Returns borrowing constant gamma
481 % boundary -- string
482 function gamm = getBoundaryBorrowing(obj, boundary)
483 switch boundary
484 case {'w','e'}
485 gamm = obj.gamm_u;
486 case {'s','n'}
487 gamm = obj.gamm_v;
488 otherwise
489 error('No such boundary: boundary = %s',boundary);
490 end
491 end
492
493 function N = size(obj)
494 N = prod(obj.m);
495 end
496 end
497 end