Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 922:1d91c2a8aada feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor LaplCurv.interfaceNonConf. Add method for interpolation operator generation in LaplCurv
author Martin Almquist <malmquist@stanford.edu>
date Sun, 02 Dec 2018 17:10:07 -0800
parents 19c05eefc8c6
children 8f8f5ff23ead
comparison
equal deleted inserted replaced
921:19c05eefc8c6 922:1d91c2a8aada
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 if isempty(opts)
282 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary); 282 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary);
283 else 283 else
284 assertType(opts, 'struct'); 284 assertType(opts, 'struct');
285 if isfield(opts, 'I_local2neighbor') 285 if isfield(opts, 'interpolation')
286 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts); 286 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
287 else 287 else
288 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts); 288 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
289 end 289 end
290 end 290 end
327 327
328 % TODO: Make this work for curvilinear grids 328 % TODO: Make this work for curvilinear grids
329 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); 329 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
330 330
331 default_field(opts, 'tuning', 1.2); 331 default_field(opts, 'tuning', 1.2);
332 default_field(opts, 'interpolation', 'OP');
332 tuning = opts.tuning; 333 tuning = opts.tuning;
333 334
334 % u denotes the solution in the own domain 335 % u denotes the solution in the own domain
335 % v denotes the solution in the neighbour domain 336 % v denotes the solution in the neighbour domain
336 I_u2v_good = opts.I_local2neighbor.good; 337 [e_u, d_u, gamm_u, H_b_u, I_u, ~, X_u] = obj.get_boundary_ops(boundary);
337 I_u2v_bad = opts.I_local2neighbor.bad; 338 [e_v, d_v, gamm_v, H_b_v, I_v, ~, X_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
338 I_v2u_good = opts.I_neighbor2local.good; 339
339 I_v2u_bad = opts.I_neighbor2local.bad; 340 % Extract the coordinate that varies along the interface
340 341 switch boundary
341 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); 342 case {'e','w'}
342 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 343 x_u = X_u(:,2);
344 case {'s', 'n'}
345 x_u = X_u(:,1);
346 end
347
348 switch neighbour_boundary
349 case {'e','w'}
350 x_v = X_v(:,2);
351 case {'s', 'n'}
352 x_v = X_v(:,1);
353 end
354
343 Hi = obj.Hi; 355 Hi = obj.Hi;
344 a = obj.a; 356 a = obj.a;
345 357
346 u = obj; 358 u = obj;
347 v = neighbour_scheme; 359 v = neighbour_scheme;
356 368
357 tau_u = tuning * spdiag(tau_u); 369 tau_u = tuning * spdiag(tau_u);
358 tau_v = tuning * spdiag(tau_v); 370 tau_v = tuning * spdiag(tau_v);
359 beta_u = tau_v; 371 beta_u = tau_v;
360 372
373 % Build interpolation operators
374 [I_u2v_good, I_u2v_bad, I_v2u_good, I_v2u_bad] = ...
375 obj.InterpolationOperators(x_u, x_v, obj.order, opts.interpolation);
376
361 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... 377 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
362 a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ... 378 a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ...
363 a*1/2*Hi*d_u*H_b_u*e_u' + ... 379 a*1/2*Hi*d_u*H_b_u*e_u' + ...
364 -a*1/2*Hi*e_u*H_b_u*d_u'; 380 -a*1/2*Hi*e_u*H_b_u*d_u';
365 381
373 389
374 % Returns the boundary ops and sign for the boundary specified by the string boundary. 390 % Returns the boundary ops and sign for the boundary specified by the string boundary.
375 % The right boundary is considered the positive boundary 391 % The right boundary is considered the positive boundary
376 % 392 %
377 % I -- the indices of the boundary points in the grid matrix 393 % I -- the indices of the boundary points in the grid matrix
378 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) 394 % X -- coordinates along the boundary;
395 function [e, d, gamm, H_b, I, X, X_logic] = get_boundary_ops(obj, boundary)
379 396
380 % gridMatrix = zeros(obj.m(2),obj.m(1)); 397 % gridMatrix = zeros(obj.m(2),obj.m(1));
381 % gridMatrix(:) = 1:numel(gridMatrix); 398 % gridMatrix(:) = 1:numel(gridMatrix);
382 399
383 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 400 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
405 I = ind(:,end)'; 422 I = ind(:,end)';
406 otherwise 423 otherwise
407 error('No such boundary: boundary = %s',boundary); 424 error('No such boundary: boundary = %s',boundary);
408 end 425 end
409 426
427 X = obj.grid.getBoundary(boundary);
428 if isa(obg.grid, 'grid.Curvilinear'))
429 X_logic = obj.grid.logic.getBoundary(boundary);
430 else
431 % Cartesian physical coordinates are also logical coordinates
432 X_logic = X;
433 end
434
410 switch boundary 435 switch boundary
411 case {'w','e'} 436 case {'w','e'}
412 gamm = obj.gamm_u; 437 gamm = obj.gamm_u;
413 case {'s','n'} 438 case {'s','n'}
414 gamm = obj.gamm_v; 439 gamm = obj.gamm_v;
415 end 440 end
441
442 end
443
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
416 end 499 end
417 500
418 function N = size(obj) 501 function N = size(obj)
419 N = prod(obj.m); 502 N = prod(obj.m);
420 end 503 end