Mercurial > repos > public > sbplib
diff +scheme/LaplaceCurvilinear.m @ 1086:90c23fd08f59 feature/laplace_curvilinear_test
Make LaplaceCurvilinearNew the default version and remove the others
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Fri, 29 Mar 2019 14:41:36 -0700 |
parents | 1341c6cea9c1 |
children | 867307f4d80f |
line wrap: on
line diff
--- a/+scheme/LaplaceCurvilinear.m Fri Mar 29 14:24:39 2019 -0700 +++ b/+scheme/LaplaceCurvilinear.m Fri Mar 29 14:41:36 2019 -0700 @@ -2,10 +2,11 @@ properties m % Number of points in each direction, possibly a vector h % Grid spacing + dim % Number of spatial dimensions grid - order % Order accuracy for the approximation + order % Order of accuracy for the approximation a,b % Parameters of the operator @@ -22,10 +23,12 @@ % Metric coefficients J, Ji a11, a12, a22 + K x_u x_v y_u y_v + s_w, s_e, s_s, s_n % Boundary integral scale factors % Inner product and operators for logical coordinates H_u, H_v % Norms in the x and y directions @@ -36,8 +39,15 @@ du_e, dv_e du_s, dv_s du_n, dv_n + + % Borrowing constants + theta_M_u, theta_M_v + theta_R_u, theta_R_v + theta_H_u, theta_H_v + + % Temporary + lambda gamm_u, gamm_v - lambda end @@ -45,14 +55,10 @@ % Implements a*div(b*grad(u)) as a SBP scheme % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) - function obj = LaplaceCurvilinear(g ,order, a, b, opSet) + function obj = LaplaceCurvilinear(g, order, a, b, opSet) default_arg('opSet',@sbp.D2Variable); default_arg('a', 1); - default_arg('b', 1); - - if b ~=1 - error('Not implemented yet') - end + default_arg('b', @(x,y) 0*x + 1); % assert(isa(g, 'grid.Curvilinear')) if isa(a, 'function_handle') @@ -60,6 +66,17 @@ a = spdiag(a); end + if isa(b, 'function_handle') + b = grid.evalOn(g, b); + b = spdiag(b); + end + + % If b is scalar + if length(b) == 1 + b = b*speye(g.N(), g.N()); + end + + dim = 2; m = g.size(); m_u = m(1); m_v = m(2); @@ -134,29 +151,37 @@ a22 = 1./J .* (x_u.^2 + y_u.^2); lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); + K = cell(dim, dim); + K{1,1} = spdiag(y_v./J); + K{1,2} = spdiag(-y_u./J); + K{2,1} = spdiag(-x_v./J); + K{2,2} = spdiag(x_u./J); + obj.K = K; + obj.x_u = x_u; obj.x_v = x_v; obj.y_u = y_u; obj.y_v = y_v; - % Assemble full operators L_12 = spdiag(a12); - Duv = Du*L_12*Dv; - Dvu = Dv*L_12*Du; + Duv = Du*b*L_12*Dv; + Dvu = Dv*b*L_12*Du; Duu = sparse(m_tot); Dvv = sparse(m_tot); ind = grid.funcToMatrix(g, 1:m_tot); for i = 1:m_v - D = D2_u(a11(ind(:,i))); + b_a11 = b*a11; + D = D2_u(b_a11(ind(:,i))); p = ind(:,i); Duu(p,p) = D; end for i = 1:m_u - D = D2_v(a22(ind(i,:))); + b_a22 = b*a22; + D = D2_v(b_a22(ind(i,:))); p = ind(i,:); Dvv(p,p) = D; end @@ -214,14 +239,29 @@ obj.h = [h_u h_v]; obj.order = order; obj.grid = g; + obj.dim = dim; obj.a = a; obj.b = b; obj.a11 = a11; obj.a12 = a12; obj.a22 = a22; + obj.s_w = spdiag(s_w); + obj.s_e = spdiag(s_e); + obj.s_s = spdiag(s_s); + obj.s_n = spdiag(s_n); + + obj.theta_M_u = h_u*ops_u.borrowing.M.d1; + obj.theta_M_v = h_v*ops_v.borrowing.M.d1; + + obj.theta_R_u = h_u*ops_u.borrowing.R.delta_D; + obj.theta_R_v = h_v*ops_v.borrowing.R.delta_D; + + obj.theta_H_u = h_u*ops_u.borrowing.H11; + obj.theta_H_v = h_v*ops_v.borrowing.H11; + + % Temporary obj.lambda = lambda; - obj.gamm_u = h_u*ops_u.borrowing.M.d1; obj.gamm_v = h_v*ops_v.borrowing.M.d1; end @@ -238,34 +278,46 @@ default_arg('type','neumann'); default_arg('parameter', []); - [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + e = obj.getBoundaryOperator('e', boundary); + d = obj.getBoundaryOperator('d', boundary); H_b = obj.getBoundaryQuadrature(boundary); - gamm = obj.getBoundaryBorrowing(boundary); + s_b = obj.getBoundaryScaling(boundary); + [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary); + m = obj.getBoundaryNumber(boundary); + + K = obj.K; + J = obj.J; + Hi = obj.Hi; + a = obj.a; + b_b = e'*obj.b*e; + 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; + sigma = 0*b_b; + for i = 1:obj.dim + sigma = sigma + e'*J*K{i,m}*K{i,m}*e; + end + sigma = sigma/s_b; + tau = tuning*(1/th_R + obj.dim/th_H)*sigma; - tau1 = tuning * spdiag(-1./b1 - 1./b2); - tau2 = 1; + closure = a*Hi*d*b_b*H_b*e' ... + -a*Hi*e*tau*b_b*H_b*e'; - tau = (tau1*e + tau2*d)*H_b; - - closure = obj.a*obj.Hi*tau*e'; - penalty = -obj.a*obj.Hi*tau; + penalty = -a*Hi*d*b_b*H_b ... + +a*Hi*e*tau*b_b*H_b; - % Neumann boundary condition + % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn. case {'N','n','neumann'} tau1 = -1; tau2 = 0; tau = (tau1*e + tau2*d)*H_b; - closure = obj.a*obj.Hi*tau*d'; - penalty = -obj.a*obj.Hi*tau; + closure = a*Hi*tau*b_b*d'; + penalty = -a*Hi*tau*b_b; % Unknown, boundary condition @@ -280,7 +332,9 @@ % -- interpolation: type of interpolation, default 'none' function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) - defaultType.tuning = 1.2; + % error('Not implemented') + + defaultType.tuning = 1.0; defaultType.interpolation = 'none'; default_struct('type', defaultType); @@ -297,44 +351,77 @@ 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, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); - H_b_u = obj.getBoundaryQuadrature(boundary); - I_u = obj.getBoundaryIndices(boundary); - gamm_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); - 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; + b_b_u = e_u'*u.b*e_u; - 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; + b_b_v = e_v'*v.b*e_v; - sig1 = -1/2; - sig2 = 0; + %--- Penalty strength tau ------------- + sigma_u = 0*b_b_u; + sigma_v = 0*b_b_v; + 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 = dim*1/th_H_u*sigma_u; + tau_H_v = dim*1/th_H_v*sigma_v; - tau = (e_u*tau1 + tau2*d_u)*H_b_u; - sig = (sig1*e_u + sig2*d_u)*H_b_u; + tau = 1/4*tuning*(b_b_u*(tau_R_u + tau_H_u) + b_b_v*(tau_R_v + tau_H_v)); + %-------------------------------------- + + % Operators/coefficients that are only required from this side + Hi = u.Hi; + H_b = u.getBoundaryQuadrature(boundary); + a = u.a; - closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); - penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); + closure = 1/2*a*Hi*d_u*b_b_u*H_b*e_u' ... + -1/2*a*Hi*e_u*H_b*b_b_u*d_u' ... + -a*Hi*e_u*tau*H_b*e_u'; + + penalty = -1/2*a*Hi*d_u*b_b_u*H_b*e_v' ... + -1/2*a*Hi*e_u*H_b*b_b_v*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'); + warning('LaplaceCurvilinear: Non-conforming interface assumes that b is constant'); % User can request special interpolation operators by specifying type.interpOpSet default_field(type, 'interpOpSet', @sbp.InterpOpsOP); @@ -344,12 +431,14 @@ % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + 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, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_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); @@ -395,47 +484,13 @@ end % Returns the boundary operator op for the boundary specified by the string boundary. - % op -- string or a cell array of strings + % op -- string % boundary -- string - function varargout = getBoundaryOperator(obj, op, boundary) - - if ~iscell(op) - op = {op}; - end + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(op, {'e', 'd'}) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) - for i = 1:numel(op) - switch op{i} - case 'e' - switch boundary - case 'w' - e = obj.e_w; - case 'e' - e = obj.e_e; - case 's' - e = obj.e_s; - case 'n' - e = obj.e_n; - otherwise - error('No such boundary: boundary = %s',boundary); - end - varargout{i} = e; - - case 'd' - switch boundary - case 'w' - d = obj.d_w; - case 'e' - d = obj.d_e; - case 's' - d = obj.d_s; - case 'n' - d = obj.d_n; - otherwise - error('No such boundary: boundary = %s',boundary); - end - varargout{i} = d; - end - end + o = obj.([op, '_', boundary]); end % Returns square boundary quadrature matrix, of dimension @@ -443,18 +498,32 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + H_b = obj.(['H_', boundary]); + end + + % Returns square boundary quadrature scaling matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function s_b = getBoundaryScaling(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + s_b = obj.(['s_', boundary]); + end + + % Returns the coordinate number corresponding to the boundary + % + % boundary -- string + function m = getBoundaryNumber(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) switch boundary - case 'w' - H_b = obj.H_w; - case 'e' - H_b = obj.H_e; - case 's' - H_b = obj.H_s; - case 'n' - H_b = obj.H_n; - otherwise - error('No such boundary: boundary = %s',boundary); + case {'w', 'e'} + m = 1; + case {'s', 'n'} + m = 2; end end @@ -478,12 +547,16 @@ % Returns borrowing constant gamma % boundary -- string - function gamm = getBoundaryBorrowing(obj, boundary) + function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary) switch boundary case {'w','e'} - gamm = obj.gamm_u; + theta_H = obj.theta_H_u; + theta_M = obj.theta_M_u; + theta_R = obj.theta_R_u; case {'s','n'} - gamm = obj.gamm_v; + theta_H = obj.theta_H_v; + theta_M = obj.theta_M_v; + theta_R = obj.theta_R_v; otherwise error('No such boundary: boundary = %s',boundary); end