Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearVirtaMin.m @ 1137:2ff1f366e64a feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix minimum and correct borrowing in VirtaMin scheme.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 10 Jun 2019 13:27:29 +0200
parents eee71789f13b
children afd06a84b69c
comparison
equal deleted inserted replaced
1136:eee71789f13b 1137:2ff1f366e64a
35 du_w, dv_w 35 du_w, dv_w
36 du_e, dv_e 36 du_e, dv_e
37 du_s, dv_s 37 du_s, dv_s
38 du_n, dv_n 38 du_n, dv_n
39 gamm_u, gamm_v 39 gamm_u, gamm_v
40 h11_u, h11_v
40 lambda 41 lambda
42
43 % Number of boundary points in minumum check
44 bp
41 45
42 end 46 end
43 47
44 methods 48 methods
45 % Implements a*div(b*grad(u)) as a SBP scheme 49 % Implements a*div(b*grad(u)) as a SBP scheme
50 default_arg('a', 1); 54 default_arg('a', 1);
51 default_arg('b', 1); 55 default_arg('b', 1);
52 56
53 if b ~=1 57 if b ~=1
54 error('Not implemented yet') 58 error('Not implemented yet')
59 end
60
61 % Number of boundary points in minimum check
62 switch order
63 case 2
64 obj.bp = 2;
65 case 4
66 obj.bp = 4;
67 case 6
68 obj.bp = 7;
55 end 69 end
56 70
57 % assert(isa(g, 'grid.Curvilinear')) 71 % assert(isa(g, 'grid.Curvilinear'))
58 if isa(a, 'function_handle') 72 if isa(a, 'function_handle')
59 a = grid.evalOn(g, a); 73 a = grid.evalOn(g, a);
222 obj.a22 = a22; 236 obj.a22 = a22;
223 obj.lambda = lambda; 237 obj.lambda = lambda;
224 238
225 obj.gamm_u = h_u*ops_u.borrowing.M.d1; 239 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
226 obj.gamm_v = h_v*ops_v.borrowing.M.d1; 240 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
241
242 obj.h11_u = h_u*ops_u.borrowing.H11;
243 obj.h11_v = h_v*ops_v.borrowing.H11;
227 end 244 end
228 245
229 246
230 % Closure functions return the opertors applied to the own doamin to close the boundary 247 % 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. 248 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
238 default_arg('type','neumann'); 255 default_arg('type','neumann');
239 default_arg('parameter', []); 256 default_arg('parameter', []);
240 257
241 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); 258 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
242 H_b = obj.getBoundaryQuadrature(boundary); 259 H_b = obj.getBoundaryQuadrature(boundary);
243 gamm = obj.getBoundaryBorrowing(boundary); 260 [gamm, h11] = obj.getBoundaryBorrowing(boundary);
261
262 lambda = obj.lambda;
263 mx = obj.m(1);
264 my = obj.m(2);
265 bp = obj.bp;
266
267 lambda_mat = reshape(lambda, my, mx);
268 switch boundary
269 case 'w'
270 lambda_min = min(lambda_mat(:,1:bp), [], 2);
271 lambda_mat = repmat(lambda_min, 1, mx);
272 case 'e'
273 lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2);
274 lambda_mat = repmat(lambda_min, 1, mx);
275 case 's'
276 lambda_min = min(lambda_mat(1:bp,:), [], 1);
277 lambda_mat = repmat(lambda_min, my, 1);
278 case 'n'
279 lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1);
280 lambda_mat = repmat(lambda_min, my, 1);
281 end
282 lambda_min = lambda_mat(:);
244 283
245 switch type 284 switch type
246 % Dirichlet boundary condition 285 % Dirichlet boundary condition
247 case {'D','d','dirichlet'} 286 case {'D','d','dirichlet'}
248 tuning = 1.0; 287 tuning = 1.0;
249 288
250 b1 = gamm*obj.lambda./obj.a11.^2; 289 switch boundary
251 b2 = gamm*obj.lambda./obj.a22.^2; 290 case {'w', 'e'}
291 b1 = gamm*lambda_min./obj.a11.^2;
292 b2 = h11 *lambda./obj.a22.^2;
293 case {'s', 'n'}
294 b1 = h11 *lambda./obj.a11.^2;
295 b2 = gamm*lambda_min./obj.a22.^2;
296 end
252 297
253 tau1 = tuning * spdiag(-1./b1 - 1./b2); 298 tau1 = tuning * spdiag(-1./b1 - 1./b2);
254 tau2 = 1; 299 tau2 = 1;
255 300
256 tau = (tau1*e + tau2*d)*H_b; 301 tau = (tau1*e + tau2*d)*H_b;
279 % Fields: 324 % Fields:
280 % -- tuning: penalty strength, defaults to 1.2 325 % -- tuning: penalty strength, defaults to 1.2
281 % -- interpolation: type of interpolation, default 'none' 326 % -- interpolation: type of interpolation, default 'none'
282 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 327 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
283 328
284 defaultType.tuning = 1.2; 329 defaultType.tuning = 1.0;
285 defaultType.interpolation = 'none'; 330 defaultType.interpolation = 'none';
286 default_struct('type', defaultType); 331 default_struct('type', defaultType);
287 332
288 switch type.interpolation 333 switch type.interpolation
289 case {'none', ''} 334 case {'none', ''}
301 % u denotes the solution in the own domain 346 % u denotes the solution in the own domain
302 % v denotes the solution in the neighbour domain 347 % v denotes the solution in the neighbour domain
303 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); 348 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
304 H_b_u = obj.getBoundaryQuadrature(boundary); 349 H_b_u = obj.getBoundaryQuadrature(boundary);
305 I_u = obj.getBoundaryIndices(boundary); 350 I_u = obj.getBoundaryIndices(boundary);
306 gamm_u = obj.getBoundaryBorrowing(boundary); 351 [gamm_u, h11_u] = obj.getBoundaryBorrowing(boundary);
307 352
308 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); 353 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
309 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); 354 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
310 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); 355 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
311 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); 356 [gamm_v, h11_v] = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
312 357
313 u = obj; 358 u = obj;
314 v = neighbour_scheme; 359 v = neighbour_scheme;
315 360
316 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 361 % Minimum check for u
317 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; 362 lambda_u = obj.lambda;
318 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; 363 mx = obj.m(1);
319 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; 364 my = obj.m(2);
365 bp = obj.bp;
366
367 lambda_mat = reshape(lambda_u, my, mx);
368 switch boundary
369 case 'w'
370 lambda_min = min(lambda_mat(:,1:bp), [], 2);
371 lambda_mat = repmat(lambda_min, 1, mx);
372 case 'e'
373 lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2);
374 lambda_mat = repmat(lambda_min, 1, mx);
375 case 's'
376 lambda_min = min(lambda_mat(1:bp,:), [], 1);
377 lambda_mat = repmat(lambda_min, my, 1);
378 case 'n'
379 lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1);
380 lambda_mat = repmat(lambda_min, my, 1);
381 end
382 lambda_min_u = lambda_mat(:);
383
384 % Minimum check for v
385 lambda_v = v.lambda;
386 mx = v.m(1);
387 my = v.m(2);
388 bp = v.bp;
389
390 lambda_mat = reshape(lambda_v, my, mx);
391 switch boundary
392 case 'w'
393 lambda_min = min(lambda_mat(:,1:bp), [], 2);
394 lambda_mat = repmat(lambda_min, 1, mx);
395 case 'e'
396 lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2);
397 lambda_mat = repmat(lambda_min, 1, mx);
398 case 's'
399 lambda_min = min(lambda_mat(1:bp,:), [], 1);
400 lambda_mat = repmat(lambda_min, my, 1);
401 case 'n'
402 lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1);
403 lambda_mat = repmat(lambda_min, my, 1);
404 end
405 lambda_min_v = lambda_mat(:);
406
407 switch boundary
408 case {'w', 'e'}
409 b1_u = gamm_u*u.lambda_min_u(I_u)./u.a11(I_u).^2;
410 b2_u = h11_u *u.lambda(I_u)./u.a22(I_u).^2;
411 case {'s', 'n'}
412 b1_u = h11_u *u.lambda(I_u)./u.a11(I_u).^2;
413 b2_u = gamm_u*u.lambda_min(I_u)./u.a22(I_u).^2;
414 end
415
416 switch neighbour_boundary
417 case {'w', 'e'}
418 b1_v = gamm_v*v.lambda_min_v(I_v)./v.a11(I_v).^2;
419 b2_v = h11_v *v.lambda(I_v)./v.a22(I_v).^2;
420 case {'s', 'n'}
421 b1_v = h11_v *v.lambda(I_v)./v.a11(I_v).^2;
422 b2_v = gamm_v*v.lambda_min(I_v)./v.a22(I_v).^2;
423 end
320 424
321 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 425 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
322 tau1 = tuning * spdiag(tau1); 426 tau1 = tuning * spdiag(tau1);
323 tau2 = 1/2; 427 tau2 = 1/2;
324 428
477 end 581 end
478 end 582 end
479 583
480 % Returns borrowing constant gamma 584 % Returns borrowing constant gamma
481 % boundary -- string 585 % boundary -- string
482 function gamm = getBoundaryBorrowing(obj, boundary) 586 function [gamm, h11] = getBoundaryBorrowing(obj, boundary)
483 switch boundary 587 switch boundary
484 case {'w','e'} 588 case {'w','e'}
485 gamm = obj.gamm_u; 589 gamm = obj.gamm_u;
590 h11 = obj.h11_u;
486 case {'s','n'} 591 case {'s','n'}
487 gamm = obj.gamm_v; 592 gamm = obj.gamm_v;
593 h11 = obj.h11_v;
488 otherwise 594 otherwise
489 error('No such boundary: boundary = %s',boundary); 595 error('No such boundary: boundary = %s',boundary);
490 end 596 end
491 end 597 end
492 598