Mercurial > repos > public > sbplib

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