Mercurial > repos > public > sbplib
changeset 1067:9a858436f8fa feature/laplace_curvilinear_test
Implement new penalty strength for interface. Bugfix missing coeff a in Dirichlet penalty.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 22 Jan 2019 18:17:01 -0800 |
| parents | d64062bed5fb |
| children | e0ecce90f8cf |
| files | +scheme/LaplaceCurvilinearNewCorner.m |
| diffstat | 1 files changed, 67 insertions(+), 34 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/LaplaceCurvilinearNewCorner.m Tue Jan 22 12:53:13 2019 -0800 +++ b/+scheme/LaplaceCurvilinearNewCorner.m Tue Jan 22 18:17:01 2019 -0800 @@ -274,13 +274,13 @@ d = obj.getBoundaryOperator('d', boundary); H_b = obj.getBoundaryQuadrature(boundary); s_b = obj.getBoundaryScaling(boundary); - [th_H, th_M, th_R] = obj.getBoundaryBorrowing(boundary); + [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary); m = obj.getBoundaryNumber(boundary); K = obj.K; J = obj.J; Hi = obj.Hi; - a_b = e'*obj.a*e; + a = obj.a; switch type % Dirichlet boundary condition @@ -302,11 +302,11 @@ tau = tuning*(tau_R + tau_H); - closure = Hi*d*a_b*H_b*e' ... - -Hi*e*tau*H_b*e'; + closure = a*Hi*d*H_b*e' ... + -a*Hi*e*tau*H_b*e'; - penalty = -Hi*d*a_b*H_b ... - +Hi*e*tau*H_b; + penalty = -a*Hi*d*H_b ... + +a*Hi*e*tau*H_b; % Neumann boundary condition @@ -315,8 +315,8 @@ tau2 = 0; tau = (tau1*e + tau2*d)*H_b; - closure = obj.a*Hi*tau*d'; - penalty = -obj.a*Hi*tau; + closure = a*Hi*tau*d'; + penalty = -a*Hi*tau; % Unknown, boundary condition @@ -333,7 +333,7 @@ % error('Not implemented') - defaultType.tuning = 1.2; + defaultType.tuning = 1.0; defaultType.interpolation = 'none'; default_struct('type', defaultType); @@ -350,46 +350,79 @@ function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) tuning = type.tuning; + dim = obj.dim; % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - e_u = obj.getBoundaryOperator('e', boundary); - d_u = obj.getBoundaryOperator('d', boundary); - H_b_u = obj.getBoundaryQuadrature(boundary); - I_u = obj.getBoundaryIndices(boundary); - gamm_u = obj.getBoundaryBorrowing(boundary); - - e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); - d_v = neighbour_scheme.getBoundaryOperator('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); - 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; + % Boundary operators, u + e_u = u.getBoundaryOperator('e', boundary); + d_u = u.getBoundaryOperator('d', boundary); + gamm_u = u.getBoundaryBorrowing(boundary); + s_b_u = u.getBoundaryScaling(boundary); + [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary); + m_u = u.getBoundaryNumber(boundary); + + % Coefficients, u + K_u = u.K; + J_u = u.J; - tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); - tau1 = tuning * spdiag(tau1); - tau2 = 1/2; + % Boundary operators, v + e_v = v.getBoundaryOperator('e', neighbour_boundary); + d_v = v.getBoundaryOperator('d', neighbour_boundary); + gamm_v = v.getBoundaryBorrowing(neighbour_boundary); + s_b_v = v.getBoundaryScaling(neighbour_boundary); + [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary); + m_v = v.getBoundaryNumber(neighbour_boundary); + + % Coefficients, v + K_v = v.K; + J_v = v.J; - sig1 = -1/2; - sig2 = 0; + %--- Penalty strength tau ------------- + sigma_u = 0; + sigma_v = 0; + for i = 1:obj.dim + sigma_u = sigma_u + e_u'*J_u*K_u{i,m_u}*K_u{i,m_u}*e_u; + sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v; + end + sigma_u = sigma_u/s_b_u; + sigma_v = sigma_v/s_b_v; + + tau_R_u = 1/th_R_u*sigma_u; + tau_R_v = 1/th_R_v*sigma_v; + + tau_H_u = 1/th_H_u*sigma_u; + tau_H_u(1,1) = dim*tau_H_u(1,1); + tau_H_u(end,end) = dim*tau_H_u(end,end); - tau = (e_u*tau1 + tau2*d_u)*H_b_u; - sig = (sig1*e_u + sig2*d_u)*H_b_u; + tau_H_v = 1/th_H_v*sigma_v; + tau_H_v(1,1) = dim*tau_H_v(1,1); + tau_H_v(end,end) = dim*tau_H_v(end,end); + + tau = 1/4*tuning*(tau_R_u + tau_H_u + tau_R_v + tau_H_v); + %-------------------------------------- - closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); - penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); + % Operators/coefficients that are only required from this side + Hi = u.Hi; + H_b = u.getBoundaryQuadrature(boundary); + a = u.a; + + closure = 1/2*a*Hi*d_u*H_b*e_u' ... + -1/2*a*Hi*e_u*H_b*d_u' ... + -a*Hi*e_u*tau*H_b*e_u'; + + penalty = -1/2*a*Hi*d_u*H_b*e_v' ... + -1/2*a*Hi*e_u*H_b*d_v' ... + +a*Hi*e_u*tau*H_b*e_v'; end function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) % TODO: Make this work for curvilinear grids warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); + warning('LaplaceCurvilinear: Non-conforming interface uses Virtas penalty strength'); % User can request special interpolation operators by specifying type.interpOpSet default_field(type, 'interpOpSet', @sbp.InterpOpsOP);