Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNew.m @ 1063:e913cdb34dcb feature/laplace_curvilinear_test

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