Mercurial > repos > public > sbplib

changeset 1069:7a55a72729e6 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Copy-paste updates from LaplCurvNewCorner to LaplCurvNew, and then modify the corner part.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 23 Jan 2019 17:11:48 -0800
parents e0ecce90f8cf
children d5290a056049
files +scheme/LaplaceCurvilinearNew.m
diffstat 1 files changed, 117 insertions(+), 114 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/LaplaceCurvilinearNew.m	Tue Jan 22 19:18:24 2019 -0800
+++ b/+scheme/LaplaceCurvilinearNew.m	Wed Jan 23 17:11:48 2019 -0800
@@ -45,6 +45,10 @@
         theta_R_u, theta_R_v
         theta_H_u, theta_H_v
 
+        % Temporary
+        lambda
+        gamm_u, gamm_v
+
     end
 
     methods
@@ -54,11 +58,7 @@
         function obj = LaplaceCurvilinearNew(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')
@@ -66,6 +66,16 @@
                 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);
@@ -139,7 +149,7 @@
             a11 =  1./J .* (x_v.^2  + y_v.^2);
             a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
             a22 =  1./J .* (x_u.^2  + y_u.^2);
-            % lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
+            lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
 
             K = cell(dim, dim);
             K{1,1} = spdiag(y_v./J);
@@ -155,21 +165,23 @@
 
             % 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
@@ -234,7 +246,6 @@
             obj.a11 = a11;
             obj.a12 = a12;
             obj.a22 = a22;
-            % obj.lambda = lambda;
             obj.s_w = spdiag(s_w);
             obj.s_e = spdiag(s_e);
             obj.s_s = spdiag(s_s);
@@ -248,6 +259,11 @@
 
             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
 
 
@@ -262,16 +278,18 @@
             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);
             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;
+            b_b = e'*obj.b*e;
 
             switch type
                 % Dirichlet boundary condition
@@ -283,23 +301,28 @@
                         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;
+                    % tau = tuning*(1/th_R + obj.dim/th_H)*sigma;
+
+                    tau_R = 1/th_R*sigma;
+                    tau_H = obj.dim*1/th_H*sigma;
 
-                    closure = Hi*d*a_b*H_b*e' ...
-                             -Hi*e*tau*H_b*e';
+                    tau = tuning*(tau_R + tau_H);
 
-                    penalty = -Hi*d*a_b*H_b ...
-                              +Hi*e*tau*H_b;
+                    closure = a*Hi*d*b_b*H_b*e' ...
+                             -a*Hi*e*tau*b_b*H_b*e';
+
+                    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*Hi*tau*d';
-                    penalty = -obj.a*Hi*tau;
+                    closure =  a*Hi*tau*b_b*d';
+                    penalty = -a*Hi*tau*b_b;
 
 
                 % Unknown, boundary condition
@@ -314,9 +337,9 @@
         %          -- interpolation:    type of interpolation, default 'none'
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 
-            error('Not implemented')
+            % error('Not implemented')
 
-            defaultType.tuning = 1.2;
+            defaultType.tuning = 1.0;
             defaultType.interpolation = 'none';
             default_struct('type', defaultType);
 
@@ -333,44 +356,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;
+            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 = 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);
@@ -380,12 +436,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);
@@ -431,47 +489,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
@@ -479,19 +503,9 @@
         %
         % boundary -- string
         function H_b = getBoundaryQuadrature(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);
-            end
+            H_b = obj.(['H_', boundary]);
         end
 
         % Returns square boundary quadrature scaling matrix, of dimension
@@ -499,33 +513,22 @@
         %
         % boundary -- string
         function s_b = getBoundaryScaling(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
 
-            switch boundary
-                case 'w'
-                    s_b = obj.s_w;
-                case 'e'
-                    s_b = obj.s_e;
-                case 's'
-                    s_b = obj.s_s;
-                case 'n'
-                    s_b = obj.s_n;
-                otherwise
-                    error('No such boundary: boundary = %s',boundary);
-            end
+            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', 'e'}
                     m = 1;
                 case {'s', 'n'}
                     m = 2;
-                otherwise
-                    error('No such boundary: boundary = %s',boundary);
             end
         end