Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.m @ 369:9d1fc984f40d feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added some comments and cleaned up the code a little
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 26 Jan 2017 09:57:24 +0100
parents 53abf04f5e4e
children 0fd6561964b0
comparison
equal deleted inserted replaced
368:53abf04f5e4e 369:9d1fc984f40d
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 m % Number of points in each direction, possibly a vector
4 n %size of system 4 n % Size of system
5 h % Grid spacing 5 h % Grid spacing
6 x, y, z % Grid 6 x, y, z % Grid
7 X, Y, Z% Values of x and y for each grid point 7 X, Y, Z% Values of x and y for each grid point
8 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces 8 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
9 order % Order accuracy for the approximation 9 order % Order accuracy for the approximation
10 10
11 D % non-stabalized scheme operator 11 D % non-stabalized scheme operator
12 A, B, C, E 12 A, B, C, E % Symbolic coefficient matrices
13 Aevaluated,Bevaluated,Cevaluated, Eevaluated 13 Aevaluated,Bevaluated,Cevaluated, Eevaluated
14 14
15 H % Discrete norm 15 H % Discrete norm
16 % Norms in the x, y and z directions 16 Hx, Hy, Hz % Norms in the x, y and z directions
17 Hx, Hy, Hz
18 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 17 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
19 I_x,I_y, I_z, I_N 18 I_x,I_y, I_z, I_N
20 e_w, e_e, e_s, e_n, e_b, e_t 19 e_w, e_e, e_s, e_n, e_b, e_t
21 params %parameters for the coeficient matrice 20 params % Parameters for the coeficient matrice
22 end 21 end
23 22
24 23
25 methods 24 methods
25 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu
26 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) 26 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator)
27 default_arg('E', []) 27 default_arg('E', [])
28 default_arg('operatpr',[]) 28 default_arg('operatpr',[])
29 xlim = lim{1}; 29 xlim = lim{1};
30 ylim = lim{2}; 30 ylim = lim{2};
38 obj.B = B; 38 obj.B = B;
39 obj.C = C; 39 obj.C = C;
40 obj.E = E; 40 obj.E = E;
41 m_x = m(1); 41 m_x = m(1);
42 m_y = m(2); 42 m_y = m(2);
43 m_z=m(3); 43 m_z = m(3);
44 obj.params = params; 44 obj.params = params;
45 45
46 switch operator 46 switch operator
47 case 'upwind' 47 case 'upwind'
48 ops_x = sbp.D1Upwind(m_x,xlim,order); 48 ops_x = sbp.D1Upwind(m_x,xlim,order);
56 56
57 obj.x = ops_x.x; 57 obj.x = ops_x.x;
58 obj.y = ops_y.x; 58 obj.y = ops_y.x;
59 obj.z = ops_z.x; 59 obj.z = ops_z.x;
60 60
61 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));%% Que pasa? 61 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));
62 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); 62 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1));
63 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); 63 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z);
64 64
65 obj.Yx = kr(obj.y,ones(m_z,1)); 65 obj.Yx = kr(obj.y,ones(m_z,1));
66 obj.Zx = kr(ones(m_y,1),obj.z); 66 obj.Zx = kr(ones(m_y,1),obj.z);
67
68 obj.Xy = kr(obj.x,ones(m_z,1)); 67 obj.Xy = kr(obj.x,ones(m_z,1));
69 obj.Zy = kr(ones(m_x,1),obj.z); 68 obj.Zy = kr(ones(m_x,1),obj.z);
70
71 obj.Xz = kr(obj.x,ones(m_y,1)); 69 obj.Xz = kr(obj.x,ones(m_y,1));
72 obj.Yz = kr(ones(m_z,1),obj.y); 70 obj.Yz = kr(ones(m_z,1),obj.y);
73 71
74 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); 72 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z);
75 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); 73 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z);
83 obj.I_x = I_x; 81 obj.I_x = I_x;
84 I_y = speye(m_y); 82 I_y = speye(m_y);
85 obj.I_y = I_y; 83 obj.I_y = I_y;
86 I_z = speye(m_z); 84 I_z = speye(m_z);
87 obj.I_z = I_z; 85 obj.I_z = I_z;
88 I_N=kr(I_n,I_x,I_y,I_z); 86 I_N = kr(I_n,I_x,I_y,I_z);
89 87
90 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); 88 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z);
91 obj.Hx = ops_x.H; 89 obj.Hx = ops_x.H;
92 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); 90 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z);
93 obj.Hy = ops_y.H; 91 obj.Hy = ops_y.H;
113 111
114 Ap = (obj.Aevaluated+alphaA*I_N)/2; 112 Ap = (obj.Aevaluated+alphaA*I_N)/2;
115 Am = (obj.Aevaluated-alphaA*I_N)/2; 113 Am = (obj.Aevaluated-alphaA*I_N)/2;
116 Dpx = kr(I_n, ops_x.Dp, I_y,I_z); 114 Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
117 Dmx = kr(I_n, ops_x.Dm, I_y,I_z); 115 Dmx = kr(I_n, ops_x.Dm, I_y,I_z);
118 obj.D=-Am*Dpx; 116 obj.D = -Am*Dpx;
119 temp=Ap*Dmx; 117 temp = Ap*Dmx;
120 obj.D=obj.D-temp; 118 obj.D = obj.D-temp;
121 clear Ap Am Dpx Dmx 119 clear Ap Am Dpx Dmx
122 120
123 Bp = (obj.Bevaluated+alphaB*I_N)/2; 121 Bp = (obj.Bevaluated+alphaB*I_N)/2;
124 Bm = (obj.Bevaluated-alphaB*I_N)/2; 122 Bm = (obj.Bevaluated-alphaB*I_N)/2;
125 Dpy = kr(I_n, I_x, ops_y.Dp,I_z); 123 Dpy = kr(I_n, I_x, ops_y.Dp,I_z);
126 Dmy = kr(I_n, I_x, ops_y.Dm,I_z); 124 Dmy = kr(I_n, I_x, ops_y.Dm,I_z);
127 temp=Bm*Dpy; 125 temp = Bm*Dpy;
128 obj.D=obj.D-temp; 126 obj.D = obj.D-temp;
129 temp=Bp*Dmy; 127 temp = Bp*Dmy;
130 obj.D=obj.D-temp; 128 obj.D = obj.D-temp;
131 clear Bp Bm Dpy Dmy 129 clear Bp Bm Dpy Dmy
132 130
133 131
134 Cp = (obj.Cevaluated+alphaC*I_N)/2; 132 Cp = (obj.Cevaluated+alphaC*I_N)/2;
135 Cm = (obj.Cevaluated-alphaC*I_N)/2; 133 Cm = (obj.Cevaluated-alphaC*I_N)/2;
136 Dpz = kr(I_n, I_x, I_y,ops_z.Dp); 134 Dpz = kr(I_n, I_x, I_y,ops_z.Dp);
137 Dmz = kr(I_n, I_x, I_y,ops_z.Dm); 135 Dmz = kr(I_n, I_x, I_y,ops_z.Dm);
138 136
139 temp=Cm*Dpz; 137 temp = Cm*Dpz;
140 obj.D=obj.D-temp; 138 obj.D = obj.D-temp;
141 temp=Cp*Dmz; 139 temp = Cp*Dmz;
142 obj.D=obj.D-temp; 140 obj.D = obj.D-temp;
143 clear Cp Cm Dpz Dmz 141 clear Cp Cm Dpz Dmz
144 142 obj.D = obj.D-obj.Eevaluated;
145 obj.D=obj.D-obj.Eevaluated; 143
146 144 case 'standard'
147 otherwise
148 D1_x = kr(I_n, ops_x.D1, I_y,I_z); 145 D1_x = kr(I_n, ops_x.D1, I_y,I_z);
149 D1_y = kr(I_n, I_x, ops_y.D1,I_z); 146 D1_y = kr(I_n, I_x, ops_y.D1,I_z);
150 D1_z = kr(I_n, I_x, I_y,ops_z.D1); 147 D1_z = kr(I_n, I_x, I_y,ops_z.D1);
151 obj.D=-obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; 148 obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated;
149 otherwise
150 error('Opperator not supported');
152 end 151 end
153 end 152 end
154 153
155 % Closure functions return the opertors applied to the own doamin to close the boundary 154 % Closure functions return the opertors applied to the own doamin to close the boundary
156 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 155 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
157 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 156 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
158 % type is a string specifying the type of boundary condition if there are several. 157 % type is a string specifying the type of boundary condition if there are several.
159 % data is a function returning the data that should be applied at the boundary. 158 % data is a function returning the data that should be applied at the boundary.
160 function [closure, penalty] = boundary_condition(obj,boundary,type,L) 159 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
161 default_arg('type','char'); 160 default_arg('type','char');
162 BM=boundary_matrices(obj,boundary); 161 BM = boundary_matrices(obj,boundary);
163
164 switch type 162 switch type
165 case{'c','char'} 163 case{'c','char'}
166 [closure,penalty] = boundary_condition_char(obj,BM); 164 [closure,penalty] = boundary_condition_char(obj,BM);
167 case{'general'} 165 case{'general'}
168 [closure,penalty] = boundary_condition_general(obj,BM,boundary,L); 166 [closure,penalty] = boundary_condition_general(obj,BM,boundary,L);
183 params = obj.params; 181 params = obj.params;
184 side = max(length(X),length(Y)); 182 side = max(length(X),length(Y));
185 if isa(mat,'function_handle') 183 if isa(mat,'function_handle')
186 [rows,cols] = size(mat(params,0,0,0)); 184 [rows,cols] = size(mat(params,0,0,0));
187 matVec = mat(params,X',Y',Z'); 185 matVec = mat(params,X',Y',Z');
188 matVec=sparse(matVec); 186 matVec = sparse(matVec);
189 else 187 else
190 matVec = mat; 188 matVec = mat;
191 [rows,cols]=size(matVec); 189 [rows,cols] = size(matVec);
192 side=max(length(X),length(Y)); 190 side = max(length(X),length(Y));
193 cols=cols/side; 191 cols = cols/side;
194 end 192 end
195 ret=cell(rows,cols); 193
196 194 ret = cell(rows,cols);
197 for ii=1:rows 195 for ii = 1:rows
198 for jj=1:cols 196 for jj = 1:cols
199 ret{ii,jj}=diag(matVec(ii,(jj-1)*side+1:jj*side)); 197 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
200 end 198 end
201 end 199 end
202 ret=cell2mat(ret); 200 ret = cell2mat(ret);
203 end 201 end
204 202
205 203 function [BM] = boundary_matrices(obj,boundary)
206 function [BM]=boundary_matrices(obj,boundary) 204 params = obj.params;
207 params=obj.params;
208 205
209 switch boundary 206 switch boundary
210 case {'w','W','west'} 207 case {'w','W','west'}
211 BM.e_=obj.e_w; 208 BM.e_ = obj.e_w;
212 mat=obj.A; 209 mat = obj.A;
213 BM.boundpos='l'; 210 BM.boundpos = 'l';
214 BM.Hi=obj.Hxi; 211 BM.Hi = obj.Hxi;
215 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(1),obj.Yx,obj.Zx); 212 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(1),obj.Yx,obj.Zx);
216 BM.side=length(obj.Yx); 213 BM.side = length(obj.Yx);
217 case {'e','E','east'} 214 case {'e','E','east'}
218 BM.e_=obj.e_e; 215 BM.e_ = obj.e_e;
219 mat=obj.A; 216 mat = obj.A;
220 BM.boundpos='r'; 217 BM.boundpos = 'r';
221 BM.Hi=obj.Hxi; 218 BM.Hi = obj.Hxi;
222 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(end),obj.Yx,obj.Zx); 219 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(end),obj.Yx,obj.Zx);
223 BM.side=length(obj.Yx); 220 BM.side = length(obj.Yx);
224 case {'s','S','south'} 221 case {'s','S','south'}
225 BM.e_=obj.e_s; 222 BM.e_ = obj.e_s;
226 mat=obj.B; 223 mat = obj.B;
227 BM.boundpos='l'; 224 BM.boundpos = 'l';
228 BM.Hi=obj.Hyi; 225 BM.Hi = obj.Hyi;
229 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xy,obj.Y(1),obj.Zy); 226 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(1),obj.Zy);
230 BM.side=length(obj.Xy); 227 BM.side = length(obj.Xy);
231 case {'n','N','north'} 228 case {'n','N','north'}
232 BM.e_=obj.e_n; 229 BM.e_ = obj.e_n;
233 mat=obj.B; 230 mat = obj.B;
234 BM.boundpos='r'; 231 BM.boundpos = 'r';
235 BM.Hi=obj.Hyi; 232 BM.Hi = obj.Hyi;
236 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xy,obj.Y(end),obj.Zy); 233 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(end),obj.Zy);
237 BM.side=length(obj.Xy); 234 BM.side = length(obj.Xy);
238 case{'b','B','Bottom'} 235 case{'b','B','Bottom'}
239 BM.e_=obj.e_b; 236 BM.e_ = obj.e_b;
240 mat=obj.C; 237 mat = obj.C;
241 BM.boundpos='l'; 238 BM.boundpos = 'l';
242 BM.Hi=obj.Hzi; 239 BM.Hi = obj.Hzi;
243 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(1)); 240 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(1));
244 BM.side=length(obj.Xz); 241 BM.side = length(obj.Xz);
245 case{'t','T','Top'} 242 case{'t','T','Top'}
246 BM.e_=obj.e_t; 243 BM.e_ = obj.e_t;
247 mat=obj.C; 244 mat = obj.C;
248 BM.boundpos='r'; 245 BM.boundpos = 'r';
249 BM.Hi=obj.Hzi; 246 BM.Hi = obj.Hzi;
250 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end)); 247 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end));
251 BM.side=length(obj.Xz); 248 BM.side = length(obj.Xz);
252 end 249 end
253 250 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
254 BM.pos=signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); 251 end
255 end 252
256 253 % Characteristic bouyndary consitions
257
258 function [closure, penalty]=boundary_condition_char(obj,BM) 254 function [closure, penalty]=boundary_condition_char(obj,BM)
259 side = BM.side; 255 side = BM.side;
260 pos = BM.pos; 256 pos = BM.pos;
261 neg = BM.neg; 257 neg = BM.neg;
262 zeroval=BM.zeroval; 258 zeroval=BM.zeroval;
263 V = BM.V; 259 V = BM.V;
264 Vi = BM.Vi; 260 Vi = BM.Vi;
265 Hi=BM.Hi; 261 Hi = BM.Hi;
266 D=BM.D; 262 D = BM.D;
267 e_=BM.e_; 263 e_ = BM.e_;
268 264
269 switch BM.boundpos 265 switch BM.boundpos
270 case {'l'} 266 case {'l'}
271 tau = sparse(obj.n*side,pos); 267 tau = sparse(obj.n*side,pos);
272 Vi_plus = Vi(1:pos,:); 268 Vi_plus = Vi(1:pos,:);
280 closure = Hi*e_*V*tau*Vi_minus*e_'; 276 closure = Hi*e_*V*tau*Vi_minus*e_';
281 penalty = -Hi*e_*V*tau*Vi_minus; 277 penalty = -Hi*e_*V*tau*Vi_minus;
282 end 278 end
283 end 279 end
284 280
285 281 % General boundary condition in the form Lu=g(x)
286 function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L) 282 function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L)
287 side = BM.side; 283 side = BM.side;
288 pos = BM.pos; 284 pos = BM.pos;
289 neg = BM.neg; 285 neg = BM.neg;
290 zeroval=BM.zeroval; 286 zeroval=BM.zeroval;
291 V = BM.V; 287 V = BM.V;
292 Vi = BM.Vi; 288 Vi = BM.Vi;
293 Hi=BM.Hi; 289 Hi = BM.Hi;
294 D=BM.D; 290 D = BM.D;
295 e_=BM.e_; 291 e_ = BM.e_;
292
296 switch boundary 293 switch boundary
297 case {'w','W','west'} 294 case {'w','W','west'}
298 L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); 295 L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx);
299 case {'e','E','east'} 296 case {'e','E','east'}
300 L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx); 297 L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx);
301 case {'s','S','south'} 298 case {'s','S','south'}
302 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy); 299 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy);
303 case {'n','N','north'} 300 case {'n','N','north'}
304 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy); 301 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy);% General boundary condition in the form Lu=g(x)
305 case {'b','B','bottom'} 302 case {'b','B','bottom'}
306 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1)); 303 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1));
307 case {'t','T','top'} 304 case {'t','T','top'}
308 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); 305 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end));
309 end 306 end
332 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 329 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
333 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; 330 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
334 end 331 end
335 end 332 end
336 333
337 334 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
335 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
336 % [d+ ]
337 % D = [ d0 ]
338 % [ d-]
339 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D
338 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z) 340 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z)
339 params = obj.params; 341 params = obj.params;
340 syms xs ys zs 342 syms xs ys zs
341 [V, D] = eig(mat(params,xs,ys,zs)); 343 [V, D] = eig(mat(params,xs,ys,zs));
342 Vi=inv(V); 344 Vi=inv(V);
347 349
348 side = max(length(x),length(y)); 350 side = max(length(x),length(y));
349 Dret = zeros(obj.n,side*obj.n); 351 Dret = zeros(obj.n,side*obj.n);
350 Vret = zeros(obj.n,side*obj.n); 352 Vret = zeros(obj.n,side*obj.n);
351 Viret= zeros(obj.n,side*obj.n); 353 Viret= zeros(obj.n,side*obj.n);
354
352 for ii=1:obj.n 355 for ii=1:obj.n
353 for jj=1:obj.n 356 for jj=1:obj.n
354 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); 357 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
355 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); 358 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
356 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); 359 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));