Mercurial > repos > public > sbplib
changeset 1228:ef08adea56c4 feature/volcano
Utilize grid module in Hypsyst3d. Pass cartesian grid and replace some of the previous functionality in Hypsyst3d with properties/methods for cartesian grids.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 18 Nov 2019 10:31:58 -0800 |
| parents | 0652b34f9f27 |
| children | e1f9bedd64a9 |
| files | +scheme/Hypsyst3d.m |
| diffstat | 1 files changed, 71 insertions(+), 70 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Hypsyst3d.m Tue Nov 12 09:00:48 2019 -0800 +++ b/+scheme/Hypsyst3d.m Mon Nov 18 10:31:58 2019 -0800 @@ -1,14 +1,11 @@ classdef Hypsyst3d < scheme.Scheme properties - 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 Cartesian grid to function instead - x, y, z % Grid - X, Y, Z% Values of x and y for each grid point + grid Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces order % Order accuracy for the approximation + nEquations + D % non-stabilized scheme operator A, B, C, E % Symbolic coefficient matrices Aevaluated,Bevaluated,Cevaluated, Eevaluated @@ -24,23 +21,21 @@ 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, operator) + 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 + obj.grid = cartesianGrid; + xlim = obj.grid.lim{1}; + ylim = obj.grid.lim{2}; + zlim = obj.grid.lim{3}; obj.A = A; obj.B = B; obj.C = C; obj.E = E; - m_x = m(1); - m_y = m(2); - m_z = m(3); + m_x = obj.grid.m(1); + m_y = obj.grid.m(2); + m_z = obj.grid.m(3); obj.params = params; switch operator @@ -54,29 +49,29 @@ ops_z = sbp.D2Standard(m_z,zlim,order); end - obj.x = ops_x.x; - obj.y = ops_y.x; - obj.z = ops_z.x; - - 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); + x = obj.grid.x{1}; + y = obj.grid.x{2}; + z = obj.grid.x{3}; + points = obj.grid.points(); + X = points(:, 1); + Y = points(:, 2); + Z = points(:, 3); - 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.Yx = kr(y, ones(m_z,1)); + obj.Zx = kr(ones(m_y,1), z); + obj.Xy = kr(x, ones(m_z,1)); + obj.Zy = kr(ones(m_x,1), z); + obj.Xz = kr(x, ones(m_y,1)); + obj.Yz = kr(ones(m_z,1), y); - 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.Aevaluated = obj.evaluateCoefficientMatrix(A, X, Y, Z); + obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X, Y, Z); + obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X, Y, Z); + obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X, Y, Z); - obj.n = length(A(obj.params,0,0,0)); + obj.nEquations = length(A(obj.params,0,0,0)); - I_n = speye(obj.n); + I_n = speye(obj.nEquations); I_x = speye(m_x); obj.I_x = I_x; I_y = speye(m_y); @@ -99,15 +94,13 @@ 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.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))))); + alphaA = max(abs(eig(A(params, x(end), y(end), z(end))))); + alphaB = max(abs(eig(B(params, x(end), y(end), z(end))))); + alphaC = max(abs(eig(C(params, x(end), y(end), z(end))))); Ap = (obj.Aevaluated+alphaA*I_N)/2; Am = (obj.Aevaluated-alphaA*I_N)/2; @@ -191,7 +184,7 @@ end function N = size(obj) - N = obj.m; + N = obj.grid.m; end function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) @@ -219,6 +212,10 @@ function [BM] = boundary_matrices(obj,boundary) params = obj.params; + points = obj.grid.points(); + X = points(:, 1); + Y = points(:, 2); + Z = points(:, 3); switch boundary case {'w','W','west'} @@ -226,42 +223,42 @@ 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.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,X(1),obj.Yx,obj.Zx); BM.side = length(obj.Yx); case {'e','E','east'} 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.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,X(end),obj.Yx,obj.Zx); BM.side = length(obj.Yx); case {'s','S','south'} 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.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,Y(1),obj.Zy); BM.side = length(obj.Xy); case {'n','N','north'} 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.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,Y(end),obj.Zy); BM.side = length(obj.Xy); case{'b','B','Bottom'} 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.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,Z(1)); BM.side = length(obj.Xz); case{'t','T','Top'} 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.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,Z(end)); BM.side = length(obj.Xz); end BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); @@ -281,15 +278,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 @@ -307,41 +304,45 @@ D = BM.D; e_ = BM.e_; + x = obj.grid.x{1}; + y = obj.grid.x{2}; + z = obj.grid.x{3}; + switch boundary case {'w','W','west'} - L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); + L = obj.evaluateCoefficientMatrix(L,x(1),obj.Yx,obj.Zx); case {'e','E','east'} - L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx); + L = obj.evaluateCoefficientMatrix(L,x(end),obj.Yx,obj.Zx); case {'s','S','south'} - L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy); + L = obj.evaluateCoefficientMatrix(L,obj.Xy,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) + L = obj.evaluateCoefficientMatrix(L,obj.Xy,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)); + L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz, z(1)); case {'t','T','top'} - L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); + L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz, z(end)); end 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; @@ -365,12 +366,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));