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