Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNewCorner.m @ 1065:c2bd7f15da48 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add scheme where only corner points are multiplied by dim. Does not seem to decrease spectral radius much, but is stable.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 17 Jan 2019 18:57:54 -0800
parents
children d64062bed5fb
comparison
equal deleted inserted replaced
1064:1341c6cea9c1 1065:c2bd7f15da48
1 classdef LaplaceCurvilinearNewCorner < 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 = LaplaceCurvilinearNewCorner(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.0;
280
281 sigma = 0;
282 for i = 1:obj.dim
283 sigma = sigma + e'*J*K{i,m}*K{i,m}*e;
284 end
285 sigma = sigma/s_b;
286 % tau = tuning*(1/th_R + obj.dim/th_H)*sigma;
287
288 tau_R = 1/th_R*sigma;
289
290 tau_H = 1/th_H*sigma;
291 tau_H(1,1) = obj.dim*tau_H(1,1);
292 tau_H(end,end) = obj.dim*tau_H(end,end);
293
294 tau = tuning*(tau_R + tau_H);
295
296 closure = Hi*d*a_b*H_b*e' ...
297 -Hi*e*tau*H_b*e';
298
299 penalty = -Hi*d*a_b*H_b ...
300 +Hi*e*tau*H_b;
301
302
303 % Neumann boundary condition
304 case {'N','n','neumann'}
305 tau1 = -1;
306 tau2 = 0;
307 tau = (tau1*e + tau2*d)*H_b;
308
309 closure = obj.a*Hi*tau*d';
310 penalty = -obj.a*Hi*tau;
311
312
313 % Unknown, boundary condition
314 otherwise
315 error('No such boundary condition: type = %s',type);
316 end
317 end
318
319 % type Struct that specifies the interface coupling.
320 % Fields:
321 % -- tuning: penalty strength, defaults to 1.2
322 % -- interpolation: type of interpolation, default 'none'
323 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
324
325 error('Not implemented')
326
327 defaultType.tuning = 1.2;
328 defaultType.interpolation = 'none';
329 default_struct('type', defaultType);
330
331 switch type.interpolation
332 case {'none', ''}
333 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
334 case {'op','OP'}
335 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
336 otherwise
337 error('Unknown type of interpolation: %s ', type.interpolation);
338 end
339 end
340
341 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
342 tuning = type.tuning;
343
344 % u denotes the solution in the own domain
345 % v denotes the solution in the neighbour domain
346 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
347 H_b_u = obj.getBoundaryQuadrature(boundary);
348 I_u = obj.getBoundaryIndices(boundary);
349 gamm_u = obj.getBoundaryBorrowing(boundary);
350
351 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
352 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
353 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
354 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
355
356 u = obj;
357 v = neighbour_scheme;
358
359 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
360 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
361 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
362 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
363
364 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
365 tau1 = tuning * spdiag(tau1);
366 tau2 = 1/2;
367
368 sig1 = -1/2;
369 sig2 = 0;
370
371 tau = (e_u*tau1 + tau2*d_u)*H_b_u;
372 sig = (sig1*e_u + sig2*d_u)*H_b_u;
373
374 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
375 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
376 end
377
378 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
379
380 % TODO: Make this work for curvilinear grids
381 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
382
383 % User can request special interpolation operators by specifying type.interpOpSet
384 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
385 interpOpSet = type.interpOpSet;
386 tuning = type.tuning;
387
388
389 % u denotes the solution in the own domain
390 % v denotes the solution in the neighbour domain
391 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
392 H_b_u = obj.getBoundaryQuadrature(boundary);
393 I_u = obj.getBoundaryIndices(boundary);
394 gamm_u = obj.getBoundaryBorrowing(boundary);
395
396 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
397 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
398 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
399 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
400
401
402 % Find the number of grid points along the interface
403 m_u = size(e_u, 2);
404 m_v = size(e_v, 2);
405
406 Hi = obj.Hi;
407 a = obj.a;
408
409 u = obj;
410 v = neighbour_scheme;
411
412 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
413 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
414 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
415 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
416
417 tau_u = -1./(4*b1_u) -1./(4*b2_u);
418 tau_v = -1./(4*b1_v) -1./(4*b2_v);
419
420 tau_u = tuning * spdiag(tau_u);
421 tau_v = tuning * spdiag(tau_v);
422 beta_u = tau_v;
423
424 % Build interpolation operators
425 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
426 Iu2v = intOps.Iu2v;
427 Iv2u = intOps.Iv2u;
428
429 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
430 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
431 a*1/2*Hi*d_u*H_b_u*e_u' + ...
432 -a*1/2*Hi*e_u*H_b_u*d_u';
433
434 penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ...
435 -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ...
436 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
437 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
438
439 end
440
441 % Returns the boundary operator op for the boundary specified by the string boundary.
442 % op -- string or a cell array of strings
443 % boundary -- string
444 function varargout = getBoundaryOperator(obj, op, boundary)
445
446 if ~iscell(op)
447 op = {op};
448 end
449
450 for i = 1:numel(op)
451 switch op{i}
452 case 'e'
453 switch boundary
454 case 'w'
455 e = obj.e_w;
456 case 'e'
457 e = obj.e_e;
458 case 's'
459 e = obj.e_s;
460 case 'n'
461 e = obj.e_n;
462 otherwise
463 error('No such boundary: boundary = %s',boundary);
464 end
465 varargout{i} = e;
466
467 case 'd'
468 switch boundary
469 case 'w'
470 d = obj.d_w;
471 case 'e'
472 d = obj.d_e;
473 case 's'
474 d = obj.d_s;
475 case 'n'
476 d = obj.d_n;
477 otherwise
478 error('No such boundary: boundary = %s',boundary);
479 end
480 varargout{i} = d;
481 end
482 end
483 end
484
485 % Returns square boundary quadrature matrix, of dimension
486 % corresponding to the number of boundary points
487 %
488 % boundary -- string
489 function H_b = getBoundaryQuadrature(obj, boundary)
490
491 switch boundary
492 case 'w'
493 H_b = obj.H_w;
494 case 'e'
495 H_b = obj.H_e;
496 case 's'
497 H_b = obj.H_s;
498 case 'n'
499 H_b = obj.H_n;
500 otherwise
501 error('No such boundary: boundary = %s',boundary);
502 end
503 end
504
505 % Returns square boundary quadrature scaling matrix, of dimension
506 % corresponding to the number of boundary points
507 %
508 % boundary -- string
509 function s_b = getBoundaryScaling(obj, boundary)
510
511 switch boundary
512 case 'w'
513 s_b = obj.s_w;
514 case 'e'
515 s_b = obj.s_e;
516 case 's'
517 s_b = obj.s_s;
518 case 'n'
519 s_b = obj.s_n;
520 otherwise
521 error('No such boundary: boundary = %s',boundary);
522 end
523 end
524
525 % Returns the coordinate number corresponding to the boundary
526 %
527 % boundary -- string
528 function m = getBoundaryNumber(obj, boundary)
529
530 switch boundary
531 case {'w', 'e'}
532 m = 1;
533 case {'s', 'n'}
534 m = 2;
535 otherwise
536 error('No such boundary: boundary = %s',boundary);
537 end
538 end
539
540 % Returns the indices of the boundary points in the grid matrix
541 % boundary -- string
542 function I = getBoundaryIndices(obj, boundary)
543 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
544 switch boundary
545 case 'w'
546 I = ind(1,:);
547 case 'e'
548 I = ind(end,:);
549 case 's'
550 I = ind(:,1)';
551 case 'n'
552 I = ind(:,end)';
553 otherwise
554 error('No such boundary: boundary = %s',boundary);
555 end
556 end
557
558 % Returns borrowing constant gamma
559 % boundary -- string
560 function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary)
561 switch boundary
562 case {'w','e'}
563 theta_H = obj.theta_H_u;
564 theta_M = obj.theta_M_u;
565 theta_R = obj.theta_R_u;
566 case {'s','n'}
567 theta_H = obj.theta_H_v;
568 theta_M = obj.theta_M_v;
569 theta_R = obj.theta_R_v;
570 otherwise
571 error('No such boundary: boundary = %s',boundary);
572 end
573 end
574
575 function N = size(obj)
576 N = prod(obj.m);
577 end
578 end
579 end