Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.m @ 1248:10881b234f77 feature/volcano

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up hypsyst3d; use functionality from grid module and remove obsolete properties. Store operators for different coord dirs in cell arrays.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 20 Nov 2019 14:12:55 -0800
parents e1f9bedd64a9
children
comparison
equal deleted inserted replaced
1231:e1f9bedd64a9 1248:10881b234f77
1 classdef Hypsyst3d < scheme.Scheme 1 classdef Hypsyst3d < scheme.Scheme
2 properties 2 properties
3 grid 3 grid
4 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
5 order % Order accuracy for the approximation 4 order % Order accuracy for the approximation
6 5
7 nEquations 6 nEquations
8 7
9 D % non-stabilized scheme operator 8 D % non-stabilized scheme operator
29 obj.A = A; 28 obj.A = A;
30 obj.B = B; 29 obj.B = B;
31 obj.C = C; 30 obj.C = C;
32 obj.E = E; 31 obj.E = E;
33 obj.params = params; 32 obj.params = params;
34 dim = obj.grid.d; 33 assert(obj.grid.d == 3, 'Dimensions not correct')
35 assert(dim == 3, 'Dimensions not correct') 34 % Construct 1D operators for each coordinate dir
36 points = obj.grid.points(); 35 for i = 1:obj.grid.d
37 m = obj.grid.m; 36 ops{i} = opSet(obj.grid.m(i),obj.grid.lim{i},order);
38 for i = 1:dim 37 I{i} = speye(obj.grid.m(i));
39 ops{i} = opSet(m(i),obj.grid.lim{i},order); 38 end
40 x{i} = obj.grid.x{i}; 39 x = obj.grid.x; %Grid points in each coordinate dir
41 X{i} = points(:,i); 40 X = num2cell(obj.grid.points(),1); %Kroneckered grid values in each coordinate dir
42 end
43
44 obj.Yx = kr(x{2}, ones(m(3),1));
45 obj.Zx = kr(ones(m(2),1), x{3});
46 obj.Xy = kr(x{1}, ones(m(3),1));
47 obj.Zy = kr(ones(m(1),1), x{3});
48 obj.Xz = kr(x{1}, ones(m(2),1));
49 obj.Yz = kr(ones(m(3),1), x{2});
50 41
51 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, X{:}); 42 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, X{:});
52 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X{:}); 43 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X{:});
53 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X{:}); 44 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, X{:});
54 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X{:}); 45 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X{:});
55 46
56 obj.nEquations = length(A(obj.params,0,0,0)); 47 obj.nEquations = length(A(obj.params,0,0,0));
57
58 I_n = speye(obj.nEquations); 48 I_n = speye(obj.nEquations);
59 I_x = speye(m(1)); 49 I_N = kr(I_n,I{:});
60 I_y = speye(m(2)); 50
61 I_z = speye(m(3)); 51 obj.Hxi = kr(I_n, ops{1}.HI, I{2}, I{3});
62 I_N = kr(I_n,I_x,I_y,I_z);
63
64 obj.Hxi = kr(I_n, ops{1}.HI, I_y,I_z);
65 obj.Hx = ops{1}.H; 52 obj.Hx = ops{1}.H;
66 obj.Hyi = kr(I_n, I_x, ops{2}.HI,I_z); 53 obj.Hyi = kr(I_n, I{1}, ops{2}.HI, I{3});
67 obj.Hy = ops{2}.H; 54 obj.Hy = ops{2}.H;
68 obj.Hzi = kr(I_n, I_x,I_y, ops{3}.HI); 55 obj.Hzi = kr(I_n, I{1},I{2}, ops{3}.HI);
69 obj.Hz = ops{3}.H; 56 obj.Hz = ops{3}.H;
70 57
71 obj.e_w = kr(I_n, ops{1}.e_l, I_y,I_z); 58 obj.e_w = kr(I_n, ops{1}.e_l, I{2}, I{3});
72 obj.e_e = kr(I_n, ops{1}.e_r, I_y,I_z); 59 obj.e_e = kr(I_n, ops{1}.e_r, I{2}, I{3});
73 obj.e_s = kr(I_n, I_x, ops{2}.e_l,I_z); 60 obj.e_s = kr(I_n, I{1}, ops{2}.e_l, I{3});
74 obj.e_n = kr(I_n, I_x, ops{2}.e_r,I_z); 61 obj.e_n = kr(I_n, I{1}, ops{2}.e_r, I{3});
75 obj.e_b = kr(I_n, I_x, I_y, ops{3}.e_l); 62 obj.e_b = kr(I_n, I{1}, I{2}, ops{3}.e_l);
76 obj.e_t = kr(I_n, I_x, I_y, ops{3}.e_r); 63 obj.e_t = kr(I_n, I{1}, I{2}, ops{3}.e_r);
77 64
78 obj.order = order; 65 obj.order = order;
79 66
80 switch toString(opSet) 67 switch toString(opSet)
81 case 'sbp.D1Upwind' 68 case 'sbp.D1Upwind'
83 alphaB = max(abs(eig(B(params, x{1}(end), x{2}(end), x{3}(end))))); 70 alphaB = max(abs(eig(B(params, x{1}(end), x{2}(end), x{3}(end)))));
84 alphaC = max(abs(eig(C(params, x{1}(end), x{2}(end), x{3}(end))))); 71 alphaC = max(abs(eig(C(params, x{1}(end), x{2}(end), x{3}(end)))));
85 72
86 Ap = (obj.Aevaluated+alphaA*I_N)/2; 73 Ap = (obj.Aevaluated+alphaA*I_N)/2;
87 Am = (obj.Aevaluated-alphaA*I_N)/2; 74 Am = (obj.Aevaluated-alphaA*I_N)/2;
88 Dpx = kr(I_n, ops{1}.Dp, I_y,I_z); 75 Dpx = kr(I_n, ops{1}.Dp, I{2}, I{3});
89 Dmx = kr(I_n, ops{1}.Dm, I_y,I_z); 76 Dmx = kr(I_n, ops{1}.Dm, I{2}, I{3});
90 obj.D = -Am*Dpx; 77 obj.D = -Am*Dpx;
91 temp = Ap*Dmx; 78 temp = Ap*Dmx;
92 obj.D = obj.D-temp; 79 obj.D = obj.D-temp;
93 clear Ap Am Dpx Dmx 80 clear Ap Am Dpx Dmx
94 81
95 Bp = (obj.Bevaluated+alphaB*I_N)/2; 82 Bp = (obj.Bevaluated+alphaB*I_N)/2;
96 Bm = (obj.Bevaluated-alphaB*I_N)/2; 83 Bm = (obj.Bevaluated-alphaB*I_N)/2;
97 Dpy = kr(I_n, I_x, ops{2}.Dp,I_z); 84 Dpy = kr(I_n, I{1}, ops{2}.Dp, I{3});
98 Dmy = kr(I_n, I_x, ops{2}.Dm,I_z); 85 Dmy = kr(I_n, I{1}, ops{2}.Dm, I{3});
99 temp = Bm*Dpy; 86 temp = Bm*Dpy;
100 obj.D = obj.D-temp; 87 obj.D = obj.D-temp;
101 temp = Bp*Dmy; 88 temp = Bp*Dmy;
102 obj.D = obj.D-temp; 89 obj.D = obj.D-temp;
103 clear Bp Bm Dpy Dmy 90 clear Bp Bm Dpy Dmy
104 91
105 92
106 Cp = (obj.Cevaluated+alphaC*I_N)/2; 93 Cp = (obj.Cevaluated+alphaC*I_N)/2;
107 Cm = (obj.Cevaluated-alphaC*I_N)/2; 94 Cm = (obj.Cevaluated-alphaC*I_N)/2;
108 Dpz = kr(I_n, I_x, I_y,ops{3}.Dp); 95 Dpz = kr(I_n, I{1}, I{2}, ops{3}.Dp);
109 Dmz = kr(I_n, I_x, I_y,ops{3}.Dm); 96 Dmz = kr(I_n, I{1}, I{2}, ops{3}.Dm);
110 97
111 temp = Cm*Dpz; 98 temp = Cm*Dpz;
112 obj.D = obj.D-temp; 99 obj.D = obj.D-temp;
113 temp = Cp*Dmz; 100 temp = Cp*Dmz;
114 obj.D = obj.D-temp; 101 obj.D = obj.D-temp;
115 clear Cp Cm Dpz Dmz 102 clear Cp Cm Dpz Dmz
116 obj.D = obj.D-obj.Eevaluated; 103 obj.D = obj.D-obj.Eevaluated;
117 104
118 case {'sbp.D2Standard', 'sbp.D2Variable', 'sbp.D4Standard', 'sbp.D4Variable'} 105 case {'sbp.D2Standard', 'sbp.D2Variable', 'sbp.D4Standard', 'sbp.D4Variable'}
119 D1_x = kr(I_n, ops{1}.D1, I_y,I_z); 106 D1_x = kr(I_n, ops{1}.D1, I{2}, I{3});
120 D1_y = kr(I_n, I_x, ops{2}.D1,I_z); 107 D1_y = kr(I_n, I{1}, ops{2}.D1, I{3});
121 D1_z = kr(I_n, I_x, I_y,ops{3}.D1); 108 D1_z = kr(I_n, I{1}, I{2}, ops{3}.D1);
122 obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; 109 obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated;
123 otherwise 110 otherwise
124 error('Operator not supported'); 111 error('Operator not supported');
125 end 112 end
126 end 113 end
191 ret = cell2mat(ret); 178 ret = cell2mat(ret);
192 end 179 end
193 180
194 function [BM] = boundary_matrices(obj,boundary) 181 function [BM] = boundary_matrices(obj,boundary)
195 params = obj.params; 182 params = obj.params;
196 points = obj.grid.points(); 183 bp = num2cell(obj.grid.getBoundary(boundary),1); %Kroneckered boundary points
197 X = points(:, 1);
198 Y = points(:, 2);
199 Z = points(:, 3);
200
201 switch boundary 184 switch boundary
202 case {'w','W','west'} 185 case {'w'}
203 BM.e_ = obj.e_w; 186 BM.e_ = obj.e_w;
204 mat = obj.A; 187 mat = obj.A;
205 BM.boundpos = 'l'; 188 BM.boundpos = 'l';
206 BM.Hi = obj.Hxi; 189 BM.Hi = obj.Hxi;
207 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,X(1),obj.Yx,obj.Zx); 190 BM.side = length(bp{2});
208 BM.side = length(obj.Yx); 191 case {'e'}
209 case {'e','E','east'}
210 BM.e_ = obj.e_e; 192 BM.e_ = obj.e_e;
211 mat = obj.A; 193 mat = obj.A;
212 BM.boundpos = 'r'; 194 BM.boundpos = 'r';
213 BM.Hi = obj.Hxi; 195 BM.Hi = obj.Hxi;
214 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,X(end),obj.Yx,obj.Zx); 196 BM.side = length(bp{2});
215 BM.side = length(obj.Yx); 197 case {'s'}
216 case {'s','S','south'}
217 BM.e_ = obj.e_s; 198 BM.e_ = obj.e_s;
218 mat = obj.B; 199 mat = obj.B;
219 BM.boundpos = 'l'; 200 BM.boundpos = 'l';
220 BM.Hi = obj.Hyi; 201 BM.Hi = obj.Hyi;
221 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,Y(1),obj.Zy); 202 BM.side = length(bp{1});
222 BM.side = length(obj.Xy); 203 case {'n'}
223 case {'n','N','north'}
224 BM.e_ = obj.e_n; 204 BM.e_ = obj.e_n;
225 mat = obj.B; 205 mat = obj.B;
226 BM.boundpos = 'r'; 206 BM.boundpos = 'r';
227 BM.Hi = obj.Hyi; 207 BM.Hi = obj.Hyi;
228 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,Y(end),obj.Zy); 208 BM.side = length(bp{1});
229 BM.side = length(obj.Xy); 209 case{'d'}
230 case{'b','B','Bottom'}
231 BM.e_ = obj.e_b; 210 BM.e_ = obj.e_b;
232 mat = obj.C; 211 mat = obj.C;
233 BM.boundpos = 'l'; 212 BM.boundpos = 'l';
234 BM.Hi = obj.Hzi; 213 BM.Hi = obj.Hzi;
235 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,Z(1)); 214 BM.side = length(bp{1});
236 BM.side = length(obj.Xz); 215 case{'u'}
237 case{'t','T','Top'}
238 BM.e_ = obj.e_t; 216 BM.e_ = obj.e_t;
239 mat = obj.C; 217 mat = obj.C;
240 BM.boundpos = 'r'; 218 BM.boundpos = 'r';
241 BM.Hi = obj.Hzi; 219 BM.Hi = obj.Hzi;
242 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,Z(end)); 220 BM.side = length(bp{1});
243 BM.side = length(obj.Xz); 221 end
244 end 222 [BM.V,BM.Vi,BM.D,signVec] = obj.diagonalizeMatrix(mat,bp{:});
245 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); 223 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
246 end 224 end
247 225
248 % Characteristic boundary conditions 226 % Characteristic boundary conditions
249 function [closure, penalty]=boundary_condition_char(obj,BM) 227 function [closure, penalty]=boundary_condition_char(obj,BM)
283 Vi = BM.Vi; 261 Vi = BM.Vi;
284 Hi = BM.Hi; 262 Hi = BM.Hi;
285 D = BM.D; 263 D = BM.D;
286 e_ = BM.e_; 264 e_ = BM.e_;
287 265
288 x = obj.grid.x{1}; 266 bp = num2cell(obj.grid.getBoundary(boundary),1); %Kroneckered boundary points
289 y = obj.grid.x{2}; 267 L = obj.evaluateCoefficientMatrix(L,bp{:}); % General boundary condition in the form Lu=g(x)
290 z = obj.grid.x{3};
291
292 switch boundary
293 case {'w','W','west'}
294 L = obj.evaluateCoefficientMatrix(L,x(1),obj.Yx,obj.Zx);
295 case {'e','E','east'}
296 L = obj.evaluateCoefficientMatrix(L,x(end),obj.Yx,obj.Zx);
297 case {'s','S','south'}
298 L = obj.evaluateCoefficientMatrix(L,obj.Xy,y(1),obj.Zy);
299 case {'n','N','north'}
300 L = obj.evaluateCoefficientMatrix(L,obj.Xy,y(end),obj.Zy);% General boundary condition in the form Lu=g(x)
301 case {'b','B','bottom'}
302 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz, z(1));
303 case {'t','T','top'}
304 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz, z(end));
305 end
306 268
307 switch BM.boundpos 269 switch BM.boundpos
308 case {'l'} 270 case {'l'}
309 tau = sparse(obj.nEquations*side,pos); 271 tau = sparse(obj.nEquations*side,pos);
310 Vi_plus = Vi(1:pos,:); 272 Vi_plus = Vi(1:pos,:);
334 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign 296 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
335 % [d+ ] 297 % [d+ ]
336 % D = [ d0 ] 298 % D = [ d0 ]
337 % [ d-] 299 % [ d-]
338 % signVec is a vector specifying the number of positive, zero and negative eigenvalues of D 300 % signVec is a vector specifying the number of positive, zero and negative eigenvalues of D
339 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z) 301 function [V, Vi, D, signVec] = diagonalizeMatrix(obj, mat, x, y, z)
340 params = obj.params; 302 params = obj.params;
341 syms xs ys zs 303 syms xs ys zs
342 [V, D] = eig(mat(params,xs,ys,zs)); 304 [V, D] = eig(mat(params,xs,ys,zs));
343 Vi=inv(V); 305 Vi=inv(V);
344 xs = x; 306 xs = x;