Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.m @ 1228:ef08adea56c4 feature/volcano

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
1226:0652b34f9f27 1228:ef08adea56c4
1 classdef Hypsyst3d < scheme.Scheme 1 classdef Hypsyst3d < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 grid
4 n % Size of system
5 h % Grid spacing
6 grid %TODO: Abstract class requires a grid. Pass Cartesian grid to function instead
7 x, y, z % Grid
8 X, Y, Z% Values of x and y for each grid point
9 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces 4 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
10 order % Order accuracy for the approximation 5 order % Order accuracy for the approximation
6
7 nEquations
11 8
12 D % non-stabilized scheme operator 9 D % non-stabilized scheme operator
13 A, B, C, E % Symbolic coefficient matrices 10 A, B, C, E % Symbolic coefficient matrices
14 Aevaluated,Bevaluated,Cevaluated, Eevaluated 11 Aevaluated,Bevaluated,Cevaluated, Eevaluated
15 12
22 end 19 end
23 20
24 21
25 methods 22 methods
26 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu 23 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu
27 function obj = Hypsyst3d(m, lim, order, A, B, C, E, params, operator) 24 function obj = Hypsyst3d(cartesianGrid, order, A, B, C, E, params, operator)
25 assertType(cartesianGrid, 'grid.Cartesian');
28 default_arg('E', []) 26 default_arg('E', [])
29 xlim = lim{1}; 27 obj.grid = cartesianGrid;
30 ylim = lim{2}; 28 xlim = obj.grid.lim{1};
31 zlim = lim{3}; 29 ylim = obj.grid.lim{2};
32 30 zlim = obj.grid.lim{3};
33 if length(m) == 1
34 m = [m m m];
35 end
36 31
37 obj.A = A; 32 obj.A = A;
38 obj.B = B; 33 obj.B = B;
39 obj.C = C; 34 obj.C = C;
40 obj.E = E; 35 obj.E = E;
41 m_x = m(1); 36 m_x = obj.grid.m(1);
42 m_y = m(2); 37 m_y = obj.grid.m(2);
43 m_z = m(3); 38 m_z = obj.grid.m(3);
44 obj.params = params; 39 obj.params = params;
45 40
46 switch operator 41 switch operator
47 case 'upwind' 42 case 'upwind'
48 ops_x = sbp.D1Upwind(m_x,xlim,order); 43 ops_x = sbp.D1Upwind(m_x,xlim,order);
52 ops_x = sbp.D2Standard(m_x,xlim,order); 47 ops_x = sbp.D2Standard(m_x,xlim,order);
53 ops_y = sbp.D2Standard(m_y,ylim,order); 48 ops_y = sbp.D2Standard(m_y,ylim,order);
54 ops_z = sbp.D2Standard(m_z,zlim,order); 49 ops_z = sbp.D2Standard(m_z,zlim,order);
55 end 50 end
56 51
57 obj.x = ops_x.x; 52 x = obj.grid.x{1};
58 obj.y = ops_y.x; 53 y = obj.grid.x{2};
59 obj.z = ops_z.x; 54 z = obj.grid.x{3};
60 55 points = obj.grid.points();
61 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1)); 56 X = points(:, 1);
62 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); 57 Y = points(:, 2);
63 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); 58 Z = points(:, 3);
64 59
65 obj.Yx = kr(obj.y,ones(m_z,1)); 60 obj.Yx = kr(y, ones(m_z,1));
66 obj.Zx = kr(ones(m_y,1),obj.z); 61 obj.Zx = kr(ones(m_y,1), z);
67 obj.Xy = kr(obj.x,ones(m_z,1)); 62 obj.Xy = kr(x, ones(m_z,1));
68 obj.Zy = kr(ones(m_x,1),obj.z); 63 obj.Zy = kr(ones(m_x,1), z);
69 obj.Xz = kr(obj.x,ones(m_y,1)); 64 obj.Xz = kr(x, ones(m_y,1));
70 obj.Yz = kr(ones(m_z,1),obj.y); 65 obj.Yz = kr(ones(m_z,1), y);
71 66
72 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); 67 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, X, Y, Z);
73 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); 68 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X, Y, Z);
74 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); 69 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X, Y, Z);
75 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); 70 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X, Y, Z);
76 71
77 obj.n = length(A(obj.params,0,0,0)); 72 obj.nEquations = length(A(obj.params,0,0,0));
78 73
79 I_n = speye(obj.n); 74 I_n = speye(obj.nEquations);
80 I_x = speye(m_x); 75 I_x = speye(m_x);
81 obj.I_x = I_x; 76 obj.I_x = I_x;
82 I_y = speye(m_y); 77 I_y = speye(m_y);
83 obj.I_y = I_y; 78 obj.I_y = I_y;
84 I_z = speye(m_z); 79 I_z = speye(m_z);
97 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); 92 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z);
98 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); 93 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z);
99 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); 94 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l);
100 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); 95 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r);
101 96
102 obj.m = m;
103 obj.h = [ops_x.h ops_y.h ops_x.h];
104 obj.order = order; 97 obj.order = order;
105 98
106 switch operator 99 switch operator
107 case 'upwind' 100 case 'upwind'
108 alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end))))); 101 alphaA = max(abs(eig(A(params, x(end), y(end), z(end)))));
109 alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end))))); 102 alphaB = max(abs(eig(B(params, x(end), y(end), z(end)))));
110 alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end))))); 103 alphaC = max(abs(eig(C(params, x(end), y(end), z(end)))));
111 104
112 Ap = (obj.Aevaluated+alphaA*I_N)/2; 105 Ap = (obj.Aevaluated+alphaA*I_N)/2;
113 Am = (obj.Aevaluated-alphaA*I_N)/2; 106 Am = (obj.Aevaluated-alphaA*I_N)/2;
114 Dpx = kr(I_n, ops_x.Dp, I_y,I_z); 107 Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
115 Dmx = kr(I_n, ops_x.Dm, I_y,I_z); 108 Dmx = kr(I_n, ops_x.Dm, I_y,I_z);
189 function H = getBoundaryQuadrature(obj, boundary) 182 function H = getBoundaryQuadrature(obj, boundary)
190 error('Not implemented'); 183 error('Not implemented');
191 end 184 end
192 185
193 function N = size(obj) 186 function N = size(obj)
194 N = obj.m; 187 N = obj.grid.m;
195 end 188 end
196 189
197 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) 190 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z)
198 params = obj.params; 191 params = obj.params;
199 side = max(length(X),length(Y)); 192 side = max(length(X),length(Y));
217 ret = cell2mat(ret); 210 ret = cell2mat(ret);
218 end 211 end
219 212
220 function [BM] = boundary_matrices(obj,boundary) 213 function [BM] = boundary_matrices(obj,boundary)
221 params = obj.params; 214 params = obj.params;
215 points = obj.grid.points();
216 X = points(:, 1);
217 Y = points(:, 2);
218 Z = points(:, 3);
222 219
223 switch boundary 220 switch boundary
224 case {'w','W','west'} 221 case {'w','W','west'}
225 BM.e_ = obj.e_w; 222 BM.e_ = obj.e_w;
226 mat = obj.A; 223 mat = obj.A;
227 BM.boundpos = 'l'; 224 BM.boundpos = 'l';
228 BM.Hi = obj.Hxi; 225 BM.Hi = obj.Hxi;
229 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(1),obj.Yx,obj.Zx); 226 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,X(1),obj.Yx,obj.Zx);
230 BM.side = length(obj.Yx); 227 BM.side = length(obj.Yx);
231 case {'e','E','east'} 228 case {'e','E','east'}
232 BM.e_ = obj.e_e; 229 BM.e_ = obj.e_e;
233 mat = obj.A; 230 mat = obj.A;
234 BM.boundpos = 'r'; 231 BM.boundpos = 'r';
235 BM.Hi = obj.Hxi; 232 BM.Hi = obj.Hxi;
236 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(end),obj.Yx,obj.Zx); 233 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,X(end),obj.Yx,obj.Zx);
237 BM.side = length(obj.Yx); 234 BM.side = length(obj.Yx);
238 case {'s','S','south'} 235 case {'s','S','south'}
239 BM.e_ = obj.e_s; 236 BM.e_ = obj.e_s;
240 mat = obj.B; 237 mat = obj.B;
241 BM.boundpos = 'l'; 238 BM.boundpos = 'l';
242 BM.Hi = obj.Hyi; 239 BM.Hi = obj.Hyi;
243 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(1),obj.Zy); 240 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,Y(1),obj.Zy);
244 BM.side = length(obj.Xy); 241 BM.side = length(obj.Xy);
245 case {'n','N','north'} 242 case {'n','N','north'}
246 BM.e_ = obj.e_n; 243 BM.e_ = obj.e_n;
247 mat = obj.B; 244 mat = obj.B;
248 BM.boundpos = 'r'; 245 BM.boundpos = 'r';
249 BM.Hi = obj.Hyi; 246 BM.Hi = obj.Hyi;
250 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(end),obj.Zy); 247 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,Y(end),obj.Zy);
251 BM.side = length(obj.Xy); 248 BM.side = length(obj.Xy);
252 case{'b','B','Bottom'} 249 case{'b','B','Bottom'}
253 BM.e_ = obj.e_b; 250 BM.e_ = obj.e_b;
254 mat = obj.C; 251 mat = obj.C;
255 BM.boundpos = 'l'; 252 BM.boundpos = 'l';
256 BM.Hi = obj.Hzi; 253 BM.Hi = obj.Hzi;
257 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(1)); 254 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,Z(1));
258 BM.side = length(obj.Xz); 255 BM.side = length(obj.Xz);
259 case{'t','T','Top'} 256 case{'t','T','Top'}
260 BM.e_ = obj.e_t; 257 BM.e_ = obj.e_t;
261 mat = obj.C; 258 mat = obj.C;
262 BM.boundpos = 'r'; 259 BM.boundpos = 'r';
263 BM.Hi = obj.Hzi; 260 BM.Hi = obj.Hzi;
264 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end)); 261 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,Z(end));
265 BM.side = length(obj.Xz); 262 BM.side = length(obj.Xz);
266 end 263 end
267 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); 264 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
268 end 265 end
269 266
279 D = BM.D; 276 D = BM.D;
280 e_ = BM.e_; 277 e_ = BM.e_;
281 278
282 switch BM.boundpos 279 switch BM.boundpos
283 case {'l'} 280 case {'l'}
284 tau = sparse(obj.n*side,pos); 281 tau = sparse(obj.nEquations*side,pos);
285 Vi_plus = Vi(1:pos,:); 282 Vi_plus = Vi(1:pos,:);
286 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 283 tau(1:pos,:) = -abs(D(1:pos,1:pos));
287 closure = Hi*e_*V*tau*Vi_plus*e_'; 284 closure = Hi*e_*V*tau*Vi_plus*e_';
288 penalty = -Hi*e_*V*tau*Vi_plus; 285 penalty = -Hi*e_*V*tau*Vi_plus;
289 case {'r'} 286 case {'r'}
290 tau = sparse(obj.n*side,neg); 287 tau = sparse(obj.nEquations*side,neg);
291 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 288 tau((pos+zeroval)+1:obj.nEquations*side,:) = -abs(D((pos+zeroval)+1:obj.nEquations*side,(pos+zeroval)+1:obj.nEquations*side));
292 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 289 Vi_minus = Vi((pos+zeroval)+1:obj.nEquations*side,:);
293 closure = Hi*e_*V*tau*Vi_minus*e_'; 290 closure = Hi*e_*V*tau*Vi_minus*e_';
294 penalty = -Hi*e_*V*tau*Vi_minus; 291 penalty = -Hi*e_*V*tau*Vi_minus;
295 end 292 end
296 end 293 end
297 294
305 Vi = BM.Vi; 302 Vi = BM.Vi;
306 Hi = BM.Hi; 303 Hi = BM.Hi;
307 D = BM.D; 304 D = BM.D;
308 e_ = BM.e_; 305 e_ = BM.e_;
309 306
307 x = obj.grid.x{1};
308 y = obj.grid.x{2};
309 z = obj.grid.x{3};
310
310 switch boundary 311 switch boundary
311 case {'w','W','west'} 312 case {'w','W','west'}
312 L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); 313 L = obj.evaluateCoefficientMatrix(L,x(1),obj.Yx,obj.Zx);
313 case {'e','E','east'} 314 case {'e','E','east'}
314 L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx); 315 L = obj.evaluateCoefficientMatrix(L,x(end),obj.Yx,obj.Zx);
315 case {'s','S','south'} 316 case {'s','S','south'}
316 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy); 317 L = obj.evaluateCoefficientMatrix(L,obj.Xy,y(1),obj.Zy);
317 case {'n','N','north'} 318 case {'n','N','north'}
318 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy);% General boundary condition in the form Lu=g(x) 319 L = obj.evaluateCoefficientMatrix(L,obj.Xy,y(end),obj.Zy);% General boundary condition in the form Lu=g(x)
319 case {'b','B','bottom'} 320 case {'b','B','bottom'}
320 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1)); 321 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz, z(1));
321 case {'t','T','top'} 322 case {'t','T','top'}
322 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); 323 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz, z(end));
323 end 324 end
324 325
325 switch BM.boundpos 326 switch BM.boundpos
326 case {'l'} 327 case {'l'}
327 tau = sparse(obj.n*side,pos); 328 tau = sparse(obj.nEquations*side,pos);
328 Vi_plus = Vi(1:pos,:); 329 Vi_plus = Vi(1:pos,:);
329 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); 330 Vi_minus = Vi(pos+zeroval+1:obj.nEquations*side,:);
330 V_plus = V(:,1:pos); 331 V_plus = V(:,1:pos);
331 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 332 V_minus = V(:,(pos+zeroval)+1:obj.nEquations*side);
332 333
333 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 334 tau(1:pos,:) = -abs(D(1:pos,1:pos));
334 R = -inv(L*V_plus)*(L*V_minus); 335 R = -inv(L*V_plus)*(L*V_minus);
335 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 336 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
336 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; 337 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
337 case {'r'} 338 case {'r'}
338 tau = sparse(obj.n*side,neg); 339 tau = sparse(obj.nEquations*side,neg);
339 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 340 tau((pos+zeroval)+1:obj.nEquations*side,:) = -abs(D((pos+zeroval)+1:obj.nEquations*side,(pos+zeroval)+1:obj.nEquations*side));
340 Vi_plus = Vi(1:pos,:); 341 Vi_plus = Vi(1:pos,:);
341 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 342 Vi_minus = Vi((pos+zeroval)+1:obj.nEquations*side,:);
342 343
343 V_plus = V(:,1:pos); 344 V_plus = V(:,1:pos);
344 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 345 V_minus = V(:,(pos+zeroval)+1:obj.nEquations*side);
345 R = -inv(L*V_minus)*(L*V_plus); 346 R = -inv(L*V_minus)*(L*V_plus);
346 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 347 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
347 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; 348 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
348 end 349 end
349 end 350 end
363 ys = y; 364 ys = y;
364 zs = z; 365 zs = z;
365 366
366 367
367 side = max(length(x),length(y)); 368 side = max(length(x),length(y));
368 Dret = zeros(obj.n,side*obj.n); 369 Dret = zeros(obj.nEquations,side*obj.nEquations);
369 Vret = zeros(obj.n,side*obj.n); 370 Vret = zeros(obj.nEquations,side*obj.nEquations);
370 Viret= zeros(obj.n,side*obj.n); 371 Viret= zeros(obj.nEquations,side*obj.nEquations);
371 372
372 for ii=1:obj.n 373 for ii=1:obj.nEquations
373 for jj=1:obj.n 374 for jj=1:obj.nEquations
374 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); 375 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
375 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); 376 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
376 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); 377 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
377 end 378 end
378 end 379 end