Mercurial > repos > public > sbplib
changeset 1137:2ff1f366e64a feature/laplace_curvilinear_test
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 |
files | +scheme/LaplaceCurvilinearVirtaMin.m |
diffstat | 1 files changed, 117 insertions(+), 11 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/LaplaceCurvilinearVirtaMin.m Mon Jun 10 10:43:12 2019 +0200 +++ b/+scheme/LaplaceCurvilinearVirtaMin.m Mon Jun 10 13:27:29 2019 +0200 @@ -37,8 +37,12 @@ du_s, dv_s du_n, dv_n gamm_u, gamm_v + h11_u, h11_v lambda + % Number of boundary points in minumum check + bp + end methods @@ -54,6 +58,16 @@ error('Not implemented yet') end + % Number of boundary points in minimum check + switch order + case 2 + obj.bp = 2; + case 4 + obj.bp = 4; + case 6 + obj.bp = 7; + end + % assert(isa(g, 'grid.Curvilinear')) if isa(a, 'function_handle') a = grid.evalOn(g, a); @@ -224,6 +238,9 @@ obj.gamm_u = h_u*ops_u.borrowing.M.d1; obj.gamm_v = h_v*ops_v.borrowing.M.d1; + + obj.h11_u = h_u*ops_u.borrowing.H11; + obj.h11_v = h_v*ops_v.borrowing.H11; end @@ -240,15 +257,43 @@ [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); H_b = obj.getBoundaryQuadrature(boundary); - gamm = obj.getBoundaryBorrowing(boundary); + [gamm, h11] = obj.getBoundaryBorrowing(boundary); + + lambda = obj.lambda; + mx = obj.m(1); + my = obj.m(2); + bp = obj.bp; + + lambda_mat = reshape(lambda, my, mx); + switch boundary + case 'w' + lambda_min = min(lambda_mat(:,1:bp), [], 2); + lambda_mat = repmat(lambda_min, 1, mx); + case 'e' + lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2); + lambda_mat = repmat(lambda_min, 1, mx); + case 's' + lambda_min = min(lambda_mat(1:bp,:), [], 1); + lambda_mat = repmat(lambda_min, my, 1); + case 'n' + lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1); + lambda_mat = repmat(lambda_min, my, 1); + end + lambda_min = lambda_mat(:); switch type % Dirichlet boundary condition case {'D','d','dirichlet'} tuning = 1.0; - b1 = gamm*obj.lambda./obj.a11.^2; - b2 = gamm*obj.lambda./obj.a22.^2; + switch boundary + case {'w', 'e'} + b1 = gamm*lambda_min./obj.a11.^2; + b2 = h11 *lambda./obj.a22.^2; + case {'s', 'n'} + b1 = h11 *lambda./obj.a11.^2; + b2 = gamm*lambda_min./obj.a22.^2; + end tau1 = tuning * spdiag(-1./b1 - 1./b2); tau2 = 1; @@ -281,7 +326,7 @@ % -- interpolation: type of interpolation, default 'none' function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) - defaultType.tuning = 1.2; + defaultType.tuning = 1.0; defaultType.interpolation = 'none'; default_struct('type', defaultType); @@ -303,20 +348,79 @@ [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); H_b_u = obj.getBoundaryQuadrature(boundary); I_u = obj.getBoundaryIndices(boundary); - gamm_u = obj.getBoundaryBorrowing(boundary); + [gamm_u, h11_u] = obj.getBoundaryBorrowing(boundary); [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); - gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); + [gamm_v, h11_v] = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); u = obj; v = neighbour_scheme; - b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; - b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; - b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; - b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; + % Minimum check for u + lambda_u = obj.lambda; + mx = obj.m(1); + my = obj.m(2); + bp = obj.bp; + + lambda_mat = reshape(lambda_u, my, mx); + switch boundary + case 'w' + lambda_min = min(lambda_mat(:,1:bp), [], 2); + lambda_mat = repmat(lambda_min, 1, mx); + case 'e' + lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2); + lambda_mat = repmat(lambda_min, 1, mx); + case 's' + lambda_min = min(lambda_mat(1:bp,:), [], 1); + lambda_mat = repmat(lambda_min, my, 1); + case 'n' + lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1); + lambda_mat = repmat(lambda_min, my, 1); + end + lambda_min_u = lambda_mat(:); + + % Minimum check for v + lambda_v = v.lambda; + mx = v.m(1); + my = v.m(2); + bp = v.bp; + + lambda_mat = reshape(lambda_v, my, mx); + switch boundary + case 'w' + lambda_min = min(lambda_mat(:,1:bp), [], 2); + lambda_mat = repmat(lambda_min, 1, mx); + case 'e' + lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2); + lambda_mat = repmat(lambda_min, 1, mx); + case 's' + lambda_min = min(lambda_mat(1:bp,:), [], 1); + lambda_mat = repmat(lambda_min, my, 1); + case 'n' + lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1); + lambda_mat = repmat(lambda_min, my, 1); + end + lambda_min_v = lambda_mat(:); + + switch boundary + case {'w', 'e'} + b1_u = gamm_u*u.lambda_min_u(I_u)./u.a11(I_u).^2; + b2_u = h11_u *u.lambda(I_u)./u.a22(I_u).^2; + case {'s', 'n'} + b1_u = h11_u *u.lambda(I_u)./u.a11(I_u).^2; + b2_u = gamm_u*u.lambda_min(I_u)./u.a22(I_u).^2; + end + + switch neighbour_boundary + case {'w', 'e'} + b1_v = gamm_v*v.lambda_min_v(I_v)./v.a11(I_v).^2; + b2_v = h11_v *v.lambda(I_v)./v.a22(I_v).^2; + case {'s', 'n'} + b1_v = h11_v *v.lambda(I_v)./v.a11(I_v).^2; + b2_v = gamm_v*v.lambda_min(I_v)./v.a22(I_v).^2; + end tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); tau1 = tuning * spdiag(tau1); @@ -479,12 +583,14 @@ % Returns borrowing constant gamma % boundary -- string - function gamm = getBoundaryBorrowing(obj, boundary) + function [gamm, h11] = getBoundaryBorrowing(obj, boundary) switch boundary case {'w','e'} gamm = obj.gamm_u; + h11 = obj.h11_u; case {'s','n'} gamm = obj.gamm_v; + h11 = obj.h11_v; otherwise error('No such boundary: boundary = %s',boundary); end