Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 927:4291731570bb feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 03 Dec 2018 14:53:52 -0800
parents 8f8f5ff23ead
children 1c61d8fa9903
comparison
equal deleted inserted replaced
926:60a3bc79835a 927:4291731570bb
276 % opts Struct that specifies the interface coupling. 276 % opts Struct that specifies the interface coupling.
277 % Fields: 277 % Fields:
278 % -- tuning: penalty strength, defaults to 1.2 278 % -- tuning: penalty strength, defaults to 1.2
279 % -- interpolation: struct of interpolation operators (empty for conforming grids) 279 % -- interpolation: struct of interpolation operators (empty for conforming grids)
280 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 280 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
281 if isempty(opts) 281
282 defaultType.tuning = 1.2;
283 defaultType.interpolation = 'none';
284 default_struct('opts', defaultType);
285
286 switch opts.interpolation
287 case {'none', ''}
282 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary); 288 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary);
283 else 289 case {'op','OP'}
284 assertType(opts, 'struct'); 290 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
285 if isfield(opts, 'interpolation') 291 otherwise
286 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts); 292 error('Unknown type of interpolation: %s ', type.interpolation);
287 else
288 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
289 end
290 end 293 end
291 end 294 end
292 295
293 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 296 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
294 297
326 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 329 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
327 330
328 % TODO: Make this work for curvilinear grids 331 % TODO: Make this work for curvilinear grids
329 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); 332 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
330 333
331 default_field(opts, 'tuning', 1.2); 334 % User can request special interpolation operators by specifying type.interpOpSet
332 default_field(opts, 'interpolation', 'OP'); 335 default_field(opts, 'interpOpSet', @sbp.InterpOpsOP);
336 interpOpSet = opts.interpOpSet;
333 tuning = opts.tuning; 337 tuning = opts.tuning;
338
334 339
335 % u denotes the solution in the own domain 340 % u denotes the solution in the own domain
336 % v denotes the solution in the neighbour domain 341 % v denotes the solution in the neighbour domain
337 [e_u, d_u, gamm_u, H_b_u, I_u, ~, X_u] = obj.get_boundary_ops(boundary); 342 [e_u, d_u, gamm_u, H_b_u, I_u, ~, X_u] = obj.get_boundary_ops(boundary);
338 [e_v, d_v, gamm_v, H_b_v, I_v, ~, X_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 343 [e_v, d_v, gamm_v, H_b_v, I_v, ~, X_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
369 tau_u = tuning * spdiag(tau_u); 374 tau_u = tuning * spdiag(tau_u);
370 tau_v = tuning * spdiag(tau_v); 375 tau_v = tuning * spdiag(tau_v);
371 beta_u = tau_v; 376 beta_u = tau_v;
372 377
373 % Build interpolation operators 378 % Build interpolation operators
374 [I_u2v_good, I_u2v_bad, I_v2u_good, I_v2u_bad] = ... 379 intOps = interpOpSet(x_u, x_v, obj.order, neighbour_scheme.order);
375 obj.interpolationOperators(x_u, x_v, obj.order, opts.interpolation); 380 Iu2v = intOps.Iu2v;
381 Iv2u = intOps.Iv2u;
376 382
377 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... 383 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
378 a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ... 384 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
379 a*1/2*Hi*d_u*H_b_u*e_u' + ... 385 a*1/2*Hi*d_u*H_b_u*e_u' + ...
380 -a*1/2*Hi*e_u*H_b_u*d_u'; 386 -a*1/2*Hi*e_u*H_b_u*d_u';
381 387
382 penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ... 388 penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ...
383 -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ... 389 -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ...
384 -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ... 390 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
385 -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v'; 391 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
386
387 392
388 end 393 end
389 394
390 % Returns the boundary ops and sign for the boundary specified by the string boundary. 395 % Returns the boundary ops and sign for the boundary specified by the string boundary.
391 % The right boundary is considered the positive boundary 396 % The right boundary is considered the positive boundary
439 gamm = obj.gamm_v; 444 gamm = obj.gamm_v;
440 end 445 end
441 446
442 end 447 end
443 448
444 % x_u, x_v -- vectors of the coordinate that varies along the boundary
445 % interpOpSet -- string, e.g 'MC' or 'OP'.
446 % TODO: Allow for general x_u, x_v. Currently only works for 2:1 ratio.
447 function [I_u2v_good, I_u2v_bad, I_v2u_good, I_v2u_bad] = interpolationOperators(obj, x_u, x_v, order, interpOpSet)
448 default_arg(interpOpSet, 'OP');
449
450 m_u = length(x_u) - 1;
451 m_v = length(x_v) - 1;
452
453 if m_u == m_v
454 % Matching grids, no interpolation required (presumably)
455 error('LaplaceCurvilinear.interpolationOperators: m_u equals m_v, this interface should not need interpolation.')
456 elseif m_u/m_v == 2
457 % Block u is finer
458
459 switch interpOpSet
460 case 'MC'
461 interpOps = sbp.InterpMC(m_v+1, m_u+1, order, order);
462 I_u2v_good = interpOps.IF2C;
463 I_u2v_bad = interpOps.IF2C;
464 I_v2u_good = interpOps.IC2F;
465 I_v2u_bad = interpOps.IC2F;
466
467 case {'OP', 'AWW'}
468 interpOpsF2C = sbp.InterpAWW(m_v+1, m_u+1, order, order, 'F2C');
469 interpOpsC2F = sbp.InterpAWW(m_v+1, m_u+1, order, order, 'C2F');
470 I_u2v_good = interpOpsF2C.IF2C;
471 I_u2v_bad = interpOpsC2F.IF2C;
472 I_v2u_good = interpOpsC2F.IC2F;
473 I_v2u_bad = interpOpsF2C.IC2F;
474 end
475
476 elseif m_v/m_u == 2
477 % Block v is finer
478
479 switch interpOpSet
480 case 'MC'
481 interpOps = sbp.InterpMC(m_u+1, m_v+1, order, order);
482 I_u2v_good = interpOps.IC2F;
483 I_u2v_bad = interpOps.IC2F;
484 I_v2u_good = interpOps.IF2C;
485 I_v2u_bad = interpOps.IF2C;
486
487 case {'OP', 'AWW'}
488 interpOpsF2C = sbp.InterpAWW(m_u+1, m_v+1, order, order, 'F2C');
489 interpOpsC2F = sbp.InterpAWW(m_u+1, m_v+1, order, order, 'C2F');
490 I_u2v_good = interpOpsC2F.IC2F;
491 I_u2v_bad = interpOpsF2C.IC2F;
492 I_v2u_good = interpOpsF2C.IF2C;
493 I_v2u_bad = interpOpsC2F.IF2C;
494 end
495 else
496 error(sprintf('Interpolation operators for grid ratio %f have not yet been constructed', m_u/m_v));
497 end
498
499 end
500
501 function N = size(obj) 449 function N = size(obj)
502 N = prod(obj.m); 450 N = prod(obj.m);
503 end 451 end
504 end 452 end
505 end 453 end