Mercurial > repos > public > sbplib

changeset 1251:6424745e1b58 feature/volcano

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in latest changes from default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 20 Nov 2019 14:26:57 -0800
parents 10881b234f77 (diff) 8ec777fb473e (current diff)
children
files
diffstat 2 files changed, 124 insertions(+), 144 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Hypsyst3d.m	Wed Nov 20 22:24:06 2019 +0000
+++ b/+scheme/Hypsyst3d.m	Wed Nov 20 14:26:57 2019 -0800
@@ -1,117 +1,79 @@
 classdef Hypsyst3d < scheme.Scheme
     properties
-        m % Number of points in each direction, possibly a vector
-        n % Size of system
-        h % Grid spacing
-        x, y, z % Grid
-        X, Y, Z% Values of x and y for each grid point
-        Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
+        grid
         order % Order accuracy for the approximation
 
-        D % non-stabalized scheme operator
+        nEquations
+
+        D % non-stabilized scheme operator
         A, B, C, E % Symbolic coefficient matrices
         Aevaluated,Bevaluated,Cevaluated, Eevaluated
 
         H % Discrete norm
         Hx, Hy, Hz  % Norms in the x, y and z directions
         Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
-        I_x,I_y, I_z, I_N
         e_w, e_e, e_s, e_n, e_b, e_t
-        params % Parameters for the coeficient matrice
+        params % Parameters for the coefficient matrices
     end
 
 
     methods
         % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu
-        function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator)
+        function obj = Hypsyst3d(cartesianGrid, order, A, B, C, E, params, opSet)
+            assertType(cartesianGrid, 'grid.Cartesian');
             default_arg('E', [])
-            xlim =  lim{1};
-            ylim = lim{2};
-            zlim = lim{3};
-
-            if length(m) == 1
-                m = [m m m];
-            end
+            default_arg('opSet', @sbp.D2Standard)
+            obj.grid = cartesianGrid;
 
             obj.A = A;
             obj.B = B;
             obj.C = C;
             obj.E = E;
-            m_x = m(1);
-            m_y = m(2);
-            m_z = m(3);
             obj.params = params;
-
-            switch operator
-                case 'upwind'
-                    ops_x = sbp.D1Upwind(m_x,xlim,order);
-                    ops_y = sbp.D1Upwind(m_y,ylim,order);
-                    ops_z = sbp.D1Upwind(m_z,zlim,order);
-                otherwise
-                    ops_x = sbp.D2Standard(m_x,xlim,order);
-                    ops_y = sbp.D2Standard(m_y,ylim,order);
-                    ops_z = sbp.D2Standard(m_z,zlim,order);
+            assert(obj.grid.d == 3, 'Dimensions not correct')
+            % Construct 1D operators for each coordinate dir
+            for i = 1:obj.grid.d
+                ops{i} = opSet(obj.grid.m(i),obj.grid.lim{i},order);
+                I{i} = speye(obj.grid.m(i));
             end
-
-            obj.x = ops_x.x;
-            obj.y = ops_y.x;
-            obj.z = ops_z.x;
+            x = obj.grid.x; %Grid points in each coordinate dir
+            X = num2cell(obj.grid.points(),1); %Kroneckered grid values in each coordinate dir
 
-            obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));
-            obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1));
-            obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z);
-
-            obj.Yx = kr(obj.y,ones(m_z,1));
-            obj.Zx = kr(ones(m_y,1),obj.z);
-            obj.Xy = kr(obj.x,ones(m_z,1));
-            obj.Zy = kr(ones(m_x,1),obj.z);
-            obj.Xz = kr(obj.x,ones(m_y,1));
-            obj.Yz = kr(ones(m_z,1),obj.y);
+            obj.Aevaluated = obj.evaluateCoefficientMatrix(A, X{:});
+            obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X{:});
+            obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X{:});
+            obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X{:});
 
-            obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z);
-            obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z);
-            obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z);
-            obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z);
-
-            obj.n = length(A(obj.params,0,0,0));
-
-            I_n = speye(obj.n);
-            I_x = speye(m_x);
-            obj.I_x = I_x;
-            I_y = speye(m_y);
-            obj.I_y = I_y;
-            I_z = speye(m_z);
-            obj.I_z = I_z;
-            I_N = kr(I_n,I_x,I_y,I_z);
+            obj.nEquations = length(A(obj.params,0,0,0));
+            I_n = speye(obj.nEquations);
+            I_N = kr(I_n,I{:});
 
-            obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z);
-            obj.Hx = ops_x.H;
-            obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z);
-            obj.Hy = ops_y.H;
-            obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI);
-            obj.Hz = ops_z.H;
+            obj.Hxi = kr(I_n, ops{1}.HI, I{2}, I{3});
+            obj.Hx = ops{1}.H;
+            obj.Hyi = kr(I_n, I{1}, ops{2}.HI, I{3});
+            obj.Hy = ops{2}.H;
+            obj.Hzi = kr(I_n, I{1},I{2}, ops{3}.HI);
+            obj.Hz = ops{3}.H;
 
-            obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z);
-            obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z);
-            obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z);
-            obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z);
-            obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l);
-            obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r);
+            obj.e_w = kr(I_n, ops{1}.e_l, I{2}, I{3});
+            obj.e_e = kr(I_n, ops{1}.e_r, I{2}, I{3});
+            obj.e_s = kr(I_n, I{1}, ops{2}.e_l, I{3});
+            obj.e_n = kr(I_n, I{1}, ops{2}.e_r, I{3});
+            obj.e_b = kr(I_n, I{1}, I{2}, ops{3}.e_l);
+            obj.e_t = kr(I_n, I{1}, I{2}, ops{3}.e_r);
 
-            obj.m = m;
-            obj.h = [ops_x.h ops_y.h ops_x.h];
             obj.order = order;
 
-            switch operator
-                case 'upwind'
-                    alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end)))));
-                    alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end)))));
-                    alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end)))));
+            switch toString(opSet)
+                case 'sbp.D1Upwind'
+                    alphaA = max(abs(eig(A(params, x{1}(end), x{2}(end), x{3}(end)))));
+                    alphaB = max(abs(eig(B(params, x{1}(end), x{2}(end), x{3}(end)))));
+                    alphaC = max(abs(eig(C(params, x{1}(end), x{2}(end), x{3}(end)))));
 
                     Ap = (obj.Aevaluated+alphaA*I_N)/2;
                     Am = (obj.Aevaluated-alphaA*I_N)/2;
-                    Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
-                    Dmx = kr(I_n, ops_x.Dm, I_y,I_z);
+                    Dpx = kr(I_n, ops{1}.Dp, I{2}, I{3});
+                    Dmx = kr(I_n, ops{1}.Dm, I{2}, I{3});
                     obj.D = -Am*Dpx;
                     temp = Ap*Dmx;
                     obj.D = obj.D-temp;
@@ -119,8 +81,8 @@
 
                     Bp = (obj.Bevaluated+alphaB*I_N)/2;
                     Bm = (obj.Bevaluated-alphaB*I_N)/2;
-                    Dpy = kr(I_n, I_x, ops_y.Dp,I_z);
-                    Dmy = kr(I_n, I_x, ops_y.Dm,I_z);
+                    Dpy = kr(I_n, I{1}, ops{2}.Dp, I{3});
+                    Dmy = kr(I_n, I{1}, ops{2}.Dm, I{3});
                     temp = Bm*Dpy;
                     obj.D = obj.D-temp;
                     temp = Bp*Dmy;
@@ -130,8 +92,8 @@
 
                     Cp = (obj.Cevaluated+alphaC*I_N)/2;
                     Cm = (obj.Cevaluated-alphaC*I_N)/2;
-                    Dpz = kr(I_n, I_x, I_y,ops_z.Dp);
-                    Dmz = kr(I_n, I_x, I_y,ops_z.Dm);
+                    Dpz = kr(I_n, I{1}, I{2}, ops{3}.Dp);
+                    Dmz = kr(I_n, I{1}, I{2}, ops{3}.Dm);
 
                     temp = Cm*Dpz;
                     obj.D = obj.D-temp;
@@ -140,18 +102,18 @@
                     clear Cp Cm Dpz Dmz
                     obj.D = obj.D-obj.Eevaluated;
 
-                case 'standard'
-                    D1_x = kr(I_n, ops_x.D1, I_y,I_z);
-                    D1_y = kr(I_n, I_x, ops_y.D1,I_z);
-                    D1_z = kr(I_n, I_x, I_y,ops_z.D1);
+                case {'sbp.D2Standard', 'sbp.D2Variable', 'sbp.D4Standard', 'sbp.D4Variable'}
+                    D1_x = kr(I_n, ops{1}.D1, I{2}, I{3});
+                    D1_y = kr(I_n, I{1}, ops{2}.D1, I{3});
+                    D1_z = kr(I_n, I{1}, I{2}, ops{3}.D1);
                     obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated;
                 otherwise
-                    error('Opperator not supported');
+                    error('Operator not supported');
             end
         end
 
-        % Closure functions return the opertors applied to the own doamin to close the boundary
-        % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
+        % Closure functions return the operators applied to the own domain to close the boundary
+        % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
         %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
         %       type                is a string specifying the type of boundary condition if there are several.
         %       data                is a function returning the data that should be applied at the boundary.
@@ -172,8 +134,25 @@
             error('Not implemented');
         end
 
+        % TODO: Implement! This function should potentially replace boundary_matrices.
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op        -- string or a cell array of strings
+        % boundary  -- string
+        function varargout = getBoundaryOperator(obj, op, boundary)
+            error('Not implemented');
+        end
+
+        % TODO: Implement!
+        % Returns square boundary quadrature matrix, of dimension
+        % corresponding to the number of boundary points
+        %
+        % boundary -- string
+        function H = getBoundaryQuadrature(obj, boundary)
+            error('Not implemented');
+        end
+
         function N = size(obj)
-            N = obj.m;
+            N = obj.grid.m;
         end
 
         function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z)
@@ -201,55 +180,50 @@
 
         function [BM] = boundary_matrices(obj,boundary)
             params = obj.params;
-
+            bp = num2cell(obj.grid.getBoundary(boundary),1); %Kroneckered boundary points
             switch boundary
-                case {'w','W','west'}
+                case {'w'}
                     BM.e_ = obj.e_w;
                     mat = obj.A;
                     BM.boundpos = 'l';
                     BM.Hi = obj.Hxi;
-                    [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(1),obj.Yx,obj.Zx);
-                    BM.side = length(obj.Yx);
-                case {'e','E','east'}
+                    BM.side = length(bp{2});
+                case {'e'}
                     BM.e_ = obj.e_e;
                     mat = obj.A;
                     BM.boundpos = 'r';
                     BM.Hi = obj.Hxi;
-                    [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(end),obj.Yx,obj.Zx);
-                    BM.side = length(obj.Yx);
-                case {'s','S','south'}
+                    BM.side = length(bp{2});
+                case {'s'}
                     BM.e_ = obj.e_s;
                     mat = obj.B;
                     BM.boundpos = 'l';
                     BM.Hi = obj.Hyi;
-                    [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(1),obj.Zy);
-                    BM.side = length(obj.Xy);
-                case {'n','N','north'}
+                    BM.side = length(bp{1});
+                case {'n'}
                     BM.e_ = obj.e_n;
                     mat = obj.B;
                     BM.boundpos = 'r';
                     BM.Hi = obj.Hyi;
-                    [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(end),obj.Zy);
-                    BM.side = length(obj.Xy);
-                case{'b','B','Bottom'}
+                    BM.side = length(bp{1});
+                case{'d'}
                     BM.e_ = obj.e_b;
                     mat = obj.C;
                     BM.boundpos = 'l';
                     BM.Hi = obj.Hzi;
-                    [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(1));
-                    BM.side = length(obj.Xz);
-                case{'t','T','Top'}
+                    BM.side = length(bp{1});
+                case{'u'}
                     BM.e_ = obj.e_t;
                     mat = obj.C;
                     BM.boundpos = 'r';
                     BM.Hi = obj.Hzi;
-                    [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end));
-                    BM.side = length(obj.Xz);
+                    BM.side = length(bp{1});
             end
+            [BM.V,BM.Vi,BM.D,signVec] = obj.diagonalizeMatrix(mat,bp{:});
             BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
         end
 
-        % Characteristic bouyndary consitions
+        % Characteristic boundary conditions
         function [closure, penalty]=boundary_condition_char(obj,BM)
             side = BM.side;
             pos = BM.pos;
@@ -263,15 +237,15 @@
 
             switch BM.boundpos
                 case {'l'}
-                    tau = sparse(obj.n*side,pos);
+                    tau = sparse(obj.nEquations*side,pos);
                     Vi_plus = Vi(1:pos,:);
                     tau(1:pos,:) = -abs(D(1:pos,1:pos));
                     closure = Hi*e_*V*tau*Vi_plus*e_';
                     penalty = -Hi*e_*V*tau*Vi_plus;
                 case {'r'}
-                    tau = sparse(obj.n*side,neg);
-                    tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
-                    Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
+                    tau = sparse(obj.nEquations*side,neg);
+                    tau((pos+zeroval)+1:obj.nEquations*side,:) = -abs(D((pos+zeroval)+1:obj.nEquations*side,(pos+zeroval)+1:obj.nEquations*side));
+                    Vi_minus = Vi((pos+zeroval)+1:obj.nEquations*side,:);
                     closure = Hi*e_*V*tau*Vi_minus*e_';
                     penalty = -Hi*e_*V*tau*Vi_minus;
             end
@@ -289,41 +263,29 @@
             D = BM.D;
             e_ = BM.e_;
 
-            switch boundary
-                case {'w','W','west'}
-                    L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx);
-                case {'e','E','east'}
-                    L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx);
-                case {'s','S','south'}
-                    L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy);
-                case {'n','N','north'}
-                    L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy);% General boundary condition in the form Lu=g(x)
-                case {'b','B','bottom'}
-                    L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1));
-                case {'t','T','top'}
-                    L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end));
-            end
+            bp = num2cell(obj.grid.getBoundary(boundary),1); %Kroneckered boundary points
+            L = obj.evaluateCoefficientMatrix(L,bp{:}); % General boundary condition in the form Lu=g(x)
 
             switch BM.boundpos
                 case {'l'}
-                    tau = sparse(obj.n*side,pos);
+                    tau = sparse(obj.nEquations*side,pos);
                     Vi_plus = Vi(1:pos,:);
-                    Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
+                    Vi_minus = Vi(pos+zeroval+1:obj.nEquations*side,:);
                     V_plus = V(:,1:pos);
-                    V_minus = V(:,(pos+zeroval)+1:obj.n*side);
+                    V_minus = V(:,(pos+zeroval)+1:obj.nEquations*side);
 
                     tau(1:pos,:) = -abs(D(1:pos,1:pos));
                     R = -inv(L*V_plus)*(L*V_minus);
                     closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
                     penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
                 case {'r'}
-                    tau = sparse(obj.n*side,neg);
-                    tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
+                    tau = sparse(obj.nEquations*side,neg);
+                    tau((pos+zeroval)+1:obj.nEquations*side,:) = -abs(D((pos+zeroval)+1:obj.nEquations*side,(pos+zeroval)+1:obj.nEquations*side));
                     Vi_plus = Vi(1:pos,:);
-                    Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
+                    Vi_minus = Vi((pos+zeroval)+1:obj.nEquations*side,:);
 
                     V_plus = V(:,1:pos);
-                    V_minus = V(:,(pos+zeroval)+1:obj.n*side);
+                    V_minus = V(:,(pos+zeroval)+1:obj.nEquations*side);
                     R = -inv(L*V_minus)*(L*V_plus);
                     closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
                     penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
@@ -335,8 +297,8 @@
         %                                    [d+       ]
         %                               D =  [   d0    ]
         %                                    [       d-]
-        % signVec   is a vector specifying the number of possitive, zero and negative eigenvalues of D
-        function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z)
+        % signVec   is a vector specifying the number of positive, zero and negative eigenvalues of D
+        function [V, Vi, D, signVec] = diagonalizeMatrix(obj, mat, x, y, z)
             params = obj.params;
             syms xs ys zs
             [V, D] = eig(mat(params,xs,ys,zs));
@@ -347,12 +309,12 @@
 
 
             side = max(length(x),length(y));
-            Dret = zeros(obj.n,side*obj.n);
-            Vret = zeros(obj.n,side*obj.n);
-            Viret= zeros(obj.n,side*obj.n);
+            Dret = zeros(obj.nEquations,side*obj.nEquations);
+            Vret = zeros(obj.nEquations,side*obj.nEquations);
+            Viret= zeros(obj.nEquations,side*obj.nEquations);
 
-            for ii=1:obj.n
-                for jj=1:obj.n
+            for ii=1:obj.nEquations
+                for jj=1:obj.nEquations
                     Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
                     Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
                     Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
--- a/+scheme/Hypsyst3dCurve.m	Wed Nov 20 22:24:06 2019 +0000
+++ b/+scheme/Hypsyst3dCurve.m	Wed Nov 20 14:26:57 2019 -0800
@@ -3,6 +3,7 @@
         m % Number of points in each direction, possibly a vector
         n %size of system
         h % Grid spacing
+        grid %TODO: Abstract class requires a grid. Pass grid instead?
         X, Y, Z% Values of x and y for each grid point
         Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
 
@@ -268,6 +269,23 @@
             error('Not implemented');
         end
 
+        % TODO: Implement! This function should potentially replace boundary_matrices.
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op        -- string or a cell array of strings
+        % boundary  -- string
+        function varargout = getBoundaryOperator(obj, op, boundary)
+            error('Not implemented');
+        end
+
+        % TODO: Implement!
+        % Returns square boundary quadrature matrix, of dimension
+        % corresponding to the number of boundary points
+        %
+        % boundary -- string
+        function H = getBoundaryQuadrature(obj, boundary)
+            error('Not implemented');
+        end
+
         function N = size(obj)
             N = obj.m;
         end