Mercurial > repos > public > sbplib
changeset 1231:e1f9bedd64a9 feature/volcano
Add opSet to Hypsyst3d
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 18 Nov 2019 17:17:16 -0800 |
| parents | ef08adea56c4 |
| children | 10881b234f77 |
| files | +scheme/Hypsyst3d.m |
| diffstat | 1 files changed, 51 insertions(+), 70 deletions(-) [+] |
line wrap: on
line diff
diff -r ef08adea56c4 -r e1f9bedd64a9 +scheme/Hypsyst3d.m --- a/+scheme/Hypsyst3d.m Mon Nov 18 10:31:58 2019 -0800 +++ b/+scheme/Hypsyst3d.m Mon Nov 18 17:17:16 2019 -0800 @@ -13,7 +13,6 @@ 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 coefficient matrices end @@ -21,91 +20,73 @@ methods % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu - function obj = Hypsyst3d(cartesianGrid, 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', []) + default_arg('opSet', @sbp.D2Standard) 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 = obj.grid.m(1); - m_y = obj.grid.m(2); - m_z = obj.grid.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); + dim = obj.grid.d; + assert(dim == 3, 'Dimensions not correct') + points = obj.grid.points(); + m = obj.grid.m; + for i = 1:dim + ops{i} = opSet(m(i),obj.grid.lim{i},order); + x{i} = obj.grid.x{i}; + X{i} = points(:,i); end - 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(x{2}, ones(m(3),1)); + obj.Zx = kr(ones(m(2),1), x{3}); + obj.Xy = kr(x{1}, ones(m(3),1)); + obj.Zy = kr(ones(m(1),1), x{3}); + obj.Xz = kr(x{1}, ones(m(2),1)); + obj.Yz = kr(ones(m(3),1), x{2}); - 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, 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.Aevaluated = obj.evaluateCoefficientMatrix(A, X{:}); + obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X{:}); + obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X{:}); + obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X{:}); obj.nEquations = length(A(obj.params,0,0,0)); I_n = speye(obj.nEquations); - 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_x = speye(m(1)); + I_y = speye(m(2)); + I_z = speye(m(3)); I_N = kr(I_n,I_x,I_y,I_z); - 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_y,I_z); + obj.Hx = ops{1}.H; + obj.Hyi = kr(I_n, I_x, ops{2}.HI,I_z); + obj.Hy = ops{2}.H; + obj.Hzi = kr(I_n, I_x,I_y, 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_y,I_z); + obj.e_e = kr(I_n, ops{1}.e_r, I_y,I_z); + obj.e_s = kr(I_n, I_x, ops{2}.e_l,I_z); + obj.e_n = kr(I_n, I_x, ops{2}.e_r,I_z); + obj.e_b = kr(I_n, I_x, I_y, ops{3}.e_l); + obj.e_t = kr(I_n, I_x, I_y, ops{3}.e_r); obj.order = order; - switch operator - case 'upwind' - 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))))); + 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_y,I_z); + Dmx = kr(I_n, ops{1}.Dm, I_y,I_z); obj.D = -Am*Dpx; temp = Ap*Dmx; obj.D = obj.D-temp; @@ -113,8 +94,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_x, ops{2}.Dp,I_z); + Dmy = kr(I_n, I_x, ops{2}.Dm,I_z); temp = Bm*Dpy; obj.D = obj.D-temp; temp = Bp*Dmy; @@ -124,8 +105,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_x, I_y,ops{3}.Dp); + Dmz = kr(I_n, I_x, I_y,ops{3}.Dm); temp = Cm*Dpz; obj.D = obj.D-temp; @@ -134,13 +115,13 @@ 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_y,I_z); + D1_y = kr(I_n, I_x, ops{2}.D1,I_z); + D1_z = kr(I_n, I_x, I_y,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