Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearMin.m @ 1139:6bc93c091682 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement minimum check in new scheme.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 21 Jun 2019 16:27:49 +0200
parents eee71789f13b
children
comparison
equal deleted inserted replaced
1138:afd06a84b69c 1139:6bc93c091682
45 theta_R_u, theta_R_v 45 theta_R_u, theta_R_v
46 theta_H_u, theta_H_v 46 theta_H_u, theta_H_v
47 47
48 % Temporary, only used for nonconforming interfaces but should be removed. 48 % Temporary, only used for nonconforming interfaces but should be removed.
49 lambda 49 lambda
50
51 % Number of boundary points in minumum check
52 bp
50 end 53 end
51 54
52 methods 55 methods
53 % Implements a*div(b*grad(u)) as a SBP scheme 56 % 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?) 57 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
55 58
56 function obj = LaplaceCurvilinearMin(g, order, a, b, opSet) 59 function obj = LaplaceCurvilinearMin(g, order, a, b, opSet)
57 default_arg('opSet',@sbp.D2Variable); 60 default_arg('opSet',@sbp.D2Variable);
58 default_arg('a', 1); 61 default_arg('a', 1);
59 default_arg('b', @(x,y) 0*x + 1); 62 default_arg('b', @(x,y) 0*x + 1);
63
64 % Number of boundary points in minimum check
65 switch order
66 case 2
67 obj.bp = 2;
68 case 4
69 obj.bp = 4;
70 case 6
71 obj.bp = 7;
72 end
60 73
61 % assert(isa(g, 'grid.Curvilinear')) 74 % assert(isa(g, 'grid.Curvilinear'))
62 if isa(a, 'function_handle') 75 if isa(a, 'function_handle')
63 a = grid.evalOn(g, a); 76 a = grid.evalOn(g, a);
64 end 77 end
283 K = obj.K; 296 K = obj.K;
284 J = obj.J; 297 J = obj.J;
285 Hi = obj.Hi; 298 Hi = obj.Hi;
286 a = obj.a; 299 a = obj.a;
287 b_b = e'*obj.b*e; 300 b_b = e'*obj.b*e;
301 b = obj.b;
288 302
289 switch type 303 switch type
290 % Dirichlet boundary condition 304 % Dirichlet boundary condition
291 case {'D','d','dirichlet'} 305 case {'D','d','dirichlet'}
292 tuning = 1.0; 306 tuning = 1.0;
293 307
294 sigma = 0*b_b; 308 sigma = 0*b;
295 for i = 1:obj.dim 309 for i = 1:obj.dim
296 sigma = sigma + e'*J*K{i,m}*K{i,m}*e; 310 sigma = sigma + b*J*K{i,m}*K{i,m};
297 end 311 end
312
313 % Minimum check on sigma
314 mx = obj.m(1);
315 my = obj.m(2);
316 bp = obj.bp;
317 sigma_mat = reshape(diag(sigma), my, mx);
318 switch boundary
319 case 'w'
320 sigma_min = min(sigma_mat(:,1:bp), [], 2);
321 sigma_mat = repmat(sigma_min, 1, mx);
322 case 'e'
323 sigma_min = min(sigma_mat(:,end-bp+1:end), [], 2);
324 sigma_mat = repmat(sigma_min, 1, mx);
325 case 's'
326 sigma_min = min(sigma_mat(1:bp,:), [], 1);
327 sigma_mat = repmat(sigma_min, my, 1);
328 case 'n'
329 sigma_min = min(sigma_mat(end-bp+1:end,:), [], 1);
330 sigma_mat = repmat(sigma_min, my, 1);
331 end
332 sigma_min = sigma_mat(:);
333 sigma_min = spdiag(sigma_min);
334
335 % Window
336 sigma = e'*sigma*e;
337 sigma_min = e'*sigma_min*e;
338
298 sigma = sigma/s_b; 339 sigma = sigma/s_b;
299 tau = tuning*(1/th_R + obj.dim/th_H)*sigma; 340 sigma_min = sigma_min/s_b;
341
342 tau = tuning*(1/th_R*sigma/sigma_min*sigma + obj.dim/th_H*sigma);
300 343
301 closure = a*Hi*d*b_b*H_b*e' ... 344 closure = a*Hi*d*b_b*H_b*e' ...
302 -a*Hi*e*tau*b_b*H_b*e'; 345 -a*Hi*e*tau*H_b*e';
303 346
304 penalty = -a*Hi*d*b_b*H_b ... 347 penalty = -a*Hi*d*b_b*H_b ...
305 +a*Hi*e*tau*b_b*H_b; 348 +a*Hi*e*tau*H_b;
306 349
307 350
308 % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn. 351 % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn.
309 case {'N','n','neumann'} 352 case {'N','n','neumann'}
310 tau1 = -1; 353 tau1 = -1;
360 m_u = u.getBoundaryNumber(boundary); 403 m_u = u.getBoundaryNumber(boundary);
361 404
362 % Coefficients, u 405 % Coefficients, u
363 K_u = u.K; 406 K_u = u.K;
364 J_u = u.J; 407 J_u = u.J;
408 b_u = u.b;
365 b_b_u = e_u'*u.b*e_u; 409 b_b_u = e_u'*u.b*e_u;
366 410
367 % Boundary operators, v 411 % Boundary operators, v
368 e_v = v.getBoundaryOperator('e', neighbour_boundary); 412 e_v = v.getBoundaryOperator('e', neighbour_boundary);
369 d_v = v.getBoundaryOperator('d', neighbour_boundary); 413 d_v = v.getBoundaryOperator('d', neighbour_boundary);
379 end 423 end
380 424
381 % Coefficients, v 425 % Coefficients, v
382 K_v = v.K; 426 K_v = v.K;
383 J_v = v.J; 427 J_v = v.J;
428 b_v = v.b;
384 b_b_v = e_v'*v.b*e_v; 429 b_b_v = e_v'*v.b*e_v;
385 430
386 %--- Penalty strength tau ------------- 431 %--- Penalty strength tau -------------
387 sigma_u = 0*b_b_u; 432 sigma_u = 0*b_u;
388 sigma_v = 0*b_b_v; 433 sigma_v = 0*b_v;
389 for i = 1:obj.dim 434 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; 435 sigma_u = sigma_u + b_u*J_u*K_u{i,m_u}*K_u{i,m_u};
391 sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v; 436 sigma_v = sigma_v + b_v*J_v*K_v{i,m_v}*K_v{i,m_v};
392 end 437 end
438
439 %--- Minimum check on sigma_u ----
440 mx = u.m(1);
441 my = u.m(2);
442 bp = u.bp;
443 sigma_mat = reshape(diag(sigma_u), my, mx);
444 switch boundary
445 case 'w'
446 sigma_min = min(sigma_mat(:,1:bp), [], 2);
447 sigma_mat = repmat(sigma_min, 1, mx);
448 case 'e'
449 sigma_min = min(sigma_mat(:,end-bp+1:end), [], 2);
450 sigma_mat = repmat(sigma_min, 1, mx);
451 case 's'
452 sigma_min = min(sigma_mat(1:bp,:), [], 1);
453 sigma_mat = repmat(sigma_min, my, 1);
454 case 'n'
455 sigma_min = min(sigma_mat(end-bp+1:end,:), [], 1);
456 sigma_mat = repmat(sigma_min, my, 1);
457 end
458 sigma_min_u = sigma_mat(:);
459 sigma_min_u = spdiag(sigma_min_u);
460
461 % Window
462 sigma_u = e_u'*sigma_u*e_u;
463 sigma_min_u = e_u'*sigma_min_u*e_u;
464 % -------------------------------
465
466 %--- Minimum check on sigma_v ----
467 mx = v.m(1);
468 my = v.m(2);
469 bp = v.bp;
470 sigma_mat = reshape(diag(sigma_v), my, mx);
471 switch neighbour_boundary
472 case 'w'
473 sigma_min = min(sigma_mat(:,1:bp), [], 2);
474 sigma_mat = repmat(sigma_min, 1, mx);
475 case 'e'
476 sigma_min = min(sigma_mat(:,end-bp+1:end), [], 2);
477 sigma_mat = repmat(sigma_min, 1, mx);
478 case 's'
479 sigma_min = min(sigma_mat(1:bp,:), [], 1);
480 sigma_mat = repmat(sigma_min, my, 1);
481 case 'n'
482 sigma_min = min(sigma_mat(end-bp+1:end,:), [], 1);
483 sigma_mat = repmat(sigma_min, my, 1);
484 end
485 sigma_min_v = sigma_mat(:);
486 sigma_min_v = spdiag(sigma_min_v);
487
488 % Window
489 sigma_v = e_v'*sigma_v*e_v;
490 sigma_min_v = e_v'*sigma_min_v*e_v;
491 % -------------------------------
492
393 sigma_u = sigma_u/s_b_u; 493 sigma_u = sigma_u/s_b_u;
494 sigma_min_u = sigma_min_u/s_b_u;
495
394 sigma_v = sigma_v/s_b_v; 496 sigma_v = sigma_v/s_b_v;
395 497 sigma_min_v = sigma_min_v/s_b_v;
396 tau_R_u = 1/th_R_u*sigma_u; 498
397 tau_R_v = 1/th_R_v*sigma_v; 499 tau_R_u = 1/th_R_u*sigma_u/sigma_min_u*sigma_u;
500 tau_R_v = 1/th_R_v*sigma_v/sigma_min_v*sigma_v;
398 501
399 tau_H_u = dim*1/th_H_u*sigma_u; 502 tau_H_u = dim*1/th_H_u*sigma_u;
400 tau_H_v = dim*1/th_H_v*sigma_v; 503 tau_H_v = dim*1/th_H_v*sigma_v;
401 504
402 tau = 1/4*tuning*(b_b_u*(tau_R_u + tau_H_u) + b_b_v*(tau_R_v + tau_H_v)); 505 tau = 1/4*tuning*((tau_R_u + tau_H_u) + (tau_R_v + tau_H_v));
403 %-------------------------------------- 506 %--------------------------------------
404 507
405 % Operators/coefficients that are only required from this side 508 % Operators/coefficients that are only required from this side
406 Hi = u.Hi; 509 Hi = u.Hi;
407 H_b = u.getBoundaryQuadrature(boundary); 510 H_b = u.getBoundaryQuadrature(boundary);