Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2d.m @ 1072:6468a5f6ec79 feature/grids/LaplaceSquared

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 12 Feb 2019 17:12:42 +0100
parents 8d73fcdb07a5
children
comparison
equal deleted inserted replaced
1071:92cb03e64ca4 1072:6468a5f6ec79
4 n %size of system 4 n %size of system
5 h % Grid spacing 5 h % Grid spacing
6 x,y % Grid 6 x,y % Grid
7 X,Y % Values of x and y for each grid point 7 X,Y % Values of x and y for each grid point
8 order % Order accuracy for the approximation 8 order % Order accuracy for the approximation
9 9
10 D % non-stabalized scheme operator 10 D % non-stabalized scheme operator
11 A, B, E %Coefficient matrices 11 A, B, E %Coefficient matrices
12 12
13 H % Discrete norm 13 H % Discrete norm
14 % Norms in the x and y directions 14 % Norms in the x and y directions
15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
16 I_x,I_y, I_N 16 I_x,I_y, I_N
17 e_w, e_e, e_s, e_n 17 e_w, e_e, e_s, e_n
18 params %parameters for the coeficient matrice 18 params %parameters for the coeficient matrice
19 end 19 end
20 20
21 methods 21 methods
22 %Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu 22 %Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
23 function obj = Hypsyst2d(m, lim, order, A, B, E, params) 23 function obj = Hypsyst2d(m, lim, order, A, B, E, params)
24 default_arg('E', []) 24 default_arg('E', [])
25 xlim = lim{1}; 25 xlim = lim{1};
26 ylim = lim{2}; 26 ylim = lim{2};
27 27
28 if length(m) == 1 28 if length(m) == 1
29 m = [m m]; 29 m = [m m];
30 end 30 end
31 31
32 obj.A=A; 32 obj.A=A;
33 obj.B=B; 33 obj.B=B;
34 obj.E=E; 34 obj.E=E;
35 35
36 m_x = m(1); 36 m_x = m(1);
37 m_y = m(2); 37 m_y = m(2);
38 obj.params = params; 38 obj.params = params;
39 39
40 ops_x = sbp.D2Standard(m_x,xlim,order); 40 ops_x = sbp.D2Standard(m_x,xlim,order);
41 ops_y = sbp.D2Standard(m_y,ylim,order); 41 ops_y = sbp.D2Standard(m_y,ylim,order);
42 42
43 obj.x = ops_x.x; 43 obj.x = ops_x.x;
44 obj.y = ops_y.x; 44 obj.y = ops_y.x;
45 45
46 obj.X = kr(obj.x,ones(m_y,1)); 46 obj.X = kr(obj.x,ones(m_y,1));
47 obj.Y = kr(ones(m_x,1),obj.y); 47 obj.Y = kr(ones(m_x,1),obj.y);
48 48
49 Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y); 49 Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y);
50 Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y); 50 Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y);
51 Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y); 51 Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y);
52 52
53 obj.n = length(A(obj.params,0,0)); 53 obj.n = length(A(obj.params,0,0));
54 54
55 I_n = eye(obj.n);I_x = speye(m_x); 55 I_n = eye(obj.n);I_x = speye(m_x);
56 obj.I_x = I_x; 56 obj.I_x = I_x;
57 I_y = speye(m_y); 57 I_y = speye(m_y);
58 obj.I_y = I_y; 58 obj.I_y = I_y;
59 59
60 60
61 D1_x = kr(I_n, ops_x.D1, I_y); 61 D1_x = kr(I_n, ops_x.D1, I_y);
62 obj.Hxi = kr(I_n, ops_x.HI, I_y); 62 obj.Hxi = kr(I_n, ops_x.HI, I_y);
63 D1_y = kr(I_n, I_x, ops_y.D1); 63 D1_y = kr(I_n, I_x, ops_y.D1);
64 obj.Hyi = kr(I_n, I_x, ops_y.HI); 64 obj.Hyi = kr(I_n, I_x, ops_y.HI);
65 65
66 obj.e_w = kr(I_n, ops_x.e_l, I_y); 66 obj.e_w = kr(I_n, ops_x.e_l, I_y);
67 obj.e_e = kr(I_n, ops_x.e_r, I_y); 67 obj.e_e = kr(I_n, ops_x.e_r, I_y);
68 obj.e_s = kr(I_n, I_x, ops_y.e_l); 68 obj.e_s = kr(I_n, I_x, ops_y.e_l);
69 obj.e_n = kr(I_n, I_x, ops_y.e_r); 69 obj.e_n = kr(I_n, I_x, ops_y.e_r);
70 70
71 obj.m = m; 71 obj.m = m;
72 obj.h = [ops_x.h ops_y.h]; 72 obj.h = [ops_x.h ops_y.h];
73 obj.order = order; 73 obj.order = order;
74 74
75 obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated; 75 obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated;
76 76
77 end 77 end
78 78
79 % Closure functions return the opertors applied to the own doamin to close the boundary 79 % Closure functions return the opertors applied to the own doamin to close the boundary
80 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 80 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
81 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 81 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
82 % type is a string specifying the type of boundary condition if there are several. 82 % type is a string specifying the type of boundary condition if there are several.
83 % data is a function returning the data that should be applied at the boundary. 83 % data is a function returning the data that should be applied at the boundary.
90 [closure,penalty] = boundary_condition_general(obj,boundary,L); 90 [closure,penalty] = boundary_condition_general(obj,boundary,L);
91 otherwise 91 otherwise
92 error('No such boundary condition') 92 error('No such boundary condition')
93 end 93 end
94 end 94 end
95 95
96 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 96 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
97 error('An interface function does not exist yet'); 97 error('Not implemented');
98 end 98 end
99 99
100 function N = size(obj) 100 function N = size(obj)
101 N = obj.m; 101 N = obj.m;
102 end 102 end
103 103
104 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y) 104 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y)
105 params = obj.params; 105 params = obj.params;
106 106
107 if isa(mat,'function_handle') 107 if isa(mat,'function_handle')
108 [rows,cols] = size(mat(params,0,0)); 108 [rows,cols] = size(mat(params,0,0));
109 matVec = mat(params,X',Y'); 109 matVec = mat(params,X',Y');
110 matVec = sparse(matVec); 110 matVec = sparse(matVec);
111 side = max(length(X),length(Y)); 111 side = max(length(X),length(Y));
114 [rows,cols] = size(matVec); 114 [rows,cols] = size(matVec);
115 side = max(length(X),length(Y)); 115 side = max(length(X),length(Y));
116 cols = cols/side; 116 cols = cols/side;
117 end 117 end
118 ret = cell(rows,cols); 118 ret = cell(rows,cols);
119 119
120 for ii = 1:rows 120 for ii = 1:rows
121 for jj=1:cols 121 for jj=1:cols
122 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); 122 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
123 end 123 end
124 end 124 end
125 ret = cell2mat(ret); 125 ret = cell2mat(ret);
126 end 126 end
127 127
128 %Characteristic boundary conditions 128 %Characteristic boundary conditions
129 function [closure, penalty] = boundary_condition_char(obj,boundary) 129 function [closure, penalty] = boundary_condition_char(obj,boundary)
130 params = obj.params; 130 params = obj.params;
131 x = obj.x; 131 x = obj.x;
132 y = obj.y; 132 y = obj.y;
133 133
134 switch boundary 134 switch boundary
135 case {'w','W','west'} 135 case {'w','W','west'}
136 e_ = obj.e_w; 136 e_ = obj.e_w;
137 mat = obj.A; 137 mat = obj.A;
138 boundPos = 'l'; 138 boundPos = 'l';
162 side = max(length(x)); 162 side = max(length(x));
163 end 163 end
164 pos = signVec(1); 164 pos = signVec(1);
165 zeroval = signVec(2); 165 zeroval = signVec(2);
166 neg = signVec(3); 166 neg = signVec(3);
167 167
168 switch boundPos 168 switch boundPos
169 case {'l'} 169 case {'l'}
170 tau = sparse(obj.n*side,pos); 170 tau = sparse(obj.n*side,pos);
171 Vi_plus = Vi(1:pos,:); 171 Vi_plus = Vi(1:pos,:);
172 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 172 tau(1:pos,:) = -abs(D(1:pos,1:pos));
178 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 178 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
179 closure = Hi*e_*V*tau*Vi_minus*e_'; 179 closure = Hi*e_*V*tau*Vi_minus*e_';
180 penalty = -Hi*e_*V*tau*Vi_minus; 180 penalty = -Hi*e_*V*tau*Vi_minus;
181 end 181 end
182 end 182 end
183 183
184 % General boundary condition in the form Lu=g(x) 184 % General boundary condition in the form Lu=g(x)
185 function [closure,penalty] = boundary_condition_general(obj,boundary,L) 185 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
186 params = obj.params; 186 params = obj.params;
187 x = obj.x; 187 x = obj.x;
188 y = obj.y; 188 y = obj.y;
189 189 e_ = obj.getBoundaryOperator('e', boundary);
190
190 switch boundary 191 switch boundary
191 case {'w','W','west'} 192 case {'w','W','west'}
192 e_ = obj.e_w;
193 mat = obj.A; 193 mat = obj.A;
194 boundPos = 'l'; 194 boundPos = 'l';
195 Hi = obj.Hxi; 195 Hi = obj.Hxi;
196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y); 196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
197 L = obj.evaluateCoefficientMatrix(L,x(1),y); 197 L = obj.evaluateCoefficientMatrix(L,x(1),y);
198 side = max(length(y)); 198 side = max(length(y));
199 case {'e','E','east'} 199 case {'e','E','east'}
200 e_ = obj.e_e;
201 mat = obj.A; 200 mat = obj.A;
202 boundPos = 'r'; 201 boundPos = 'r';
203 Hi = obj.Hxi; 202 Hi = obj.Hxi;
204 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y); 203 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
205 L = obj.evaluateCoefficientMatrix(L,x(end),y); 204 L = obj.evaluateCoefficientMatrix(L,x(end),y);
206 side = max(length(y)); 205 side = max(length(y));
207 case {'s','S','south'} 206 case {'s','S','south'}
208 e_ = obj.e_s;
209 mat = obj.B; 207 mat = obj.B;
210 boundPos = 'l'; 208 boundPos = 'l';
211 Hi = obj.Hyi; 209 Hi = obj.Hyi;
212 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1)); 210 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
213 L = obj.evaluateCoefficientMatrix(L,x,y(1)); 211 L = obj.evaluateCoefficientMatrix(L,x,y(1));
214 side = max(length(x)); 212 side = max(length(x));
215 case {'n','N','north'} 213 case {'n','N','north'}
216 e_ = obj.e_n;
217 mat = obj.B; 214 mat = obj.B;
218 boundPos = 'r'; 215 boundPos = 'r';
219 Hi = obj.Hyi; 216 Hi = obj.Hyi;
220 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end)); 217 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
221 L = obj.evaluateCoefficientMatrix(L,x,y(end)); 218 L = obj.evaluateCoefficientMatrix(L,x,y(end));
222 side = max(length(x)); 219 side = max(length(x));
223 end 220 end
224 221
225 pos = signVec(1); 222 pos = signVec(1);
226 zeroval = signVec(2); 223 zeroval = signVec(2);
227 neg = signVec(3); 224 neg = signVec(3);
228 225
229 switch boundPos 226 switch boundPos
230 case {'l'} 227 case {'l'}
231 tau = sparse(obj.n*side,pos); 228 tau = sparse(obj.n*side,pos);
232 Vi_plus = Vi(1:pos,:); 229 Vi_plus = Vi(1:pos,:);
233 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); 230 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
234 V_plus = V(:,1:pos); 231 V_plus = V(:,1:pos);
235 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 232 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
236 233
237 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 234 tau(1:pos,:) = -abs(D(1:pos,1:pos));
238 R = -inv(L*V_plus)*(L*V_minus); 235 R = -inv(L*V_plus)*(L*V_minus);
239 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 236 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
240 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; 237 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
241 case {'r'} 238 case {'r'}
242 tau = sparse(obj.n*side,neg); 239 tau = sparse(obj.n*side,neg);
243 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 240 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
244 Vi_plus = Vi(1:pos,:); 241 Vi_plus = Vi(1:pos,:);
245 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 242 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
246 243
247 V_plus = V(:,1:pos); 244 V_plus = V(:,1:pos);
248 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 245 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
249 R = -inv(L*V_minus)*(L*V_plus); 246 R = -inv(L*V_minus)*(L*V_plus);
250 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 247 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
251 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; 248 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
252 end 249 end
253 end 250 end
254 251
255 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi 252 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
256 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign 253 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
257 % [d+ ] 254 % [d+ ]
258 % D = [ d0 ] 255 % D = [ d0 ]
259 % [ d-] 256 % [ d-]
260 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D 257 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D
261 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y) 258 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y)
262 params = obj.params; 259 params = obj.params;
263 syms xs ys 260 syms xs ys
264 [V, D]= eig(mat(params,xs,ys)); 261 [V, D]= eig(mat(params,xs,ys));
265 Vi = inv(V); 262 Vi = inv(V);
266 xs = x; 263 xs = x;
267 ys = y; 264 ys = y;
268 265
269 side = max(length(x),length(y)); 266 side = max(length(x),length(y));
270 Dret = zeros(obj.n,side*obj.n); 267 Dret = zeros(obj.n,side*obj.n);
271 Vret = zeros(obj.n,side*obj.n); 268 Vret = zeros(obj.n,side*obj.n);
272 Viret = zeros(obj.n,side*obj.n); 269 Viret = zeros(obj.n,side*obj.n);
273 270
274 for ii = 1:obj.n 271 for ii = 1:obj.n
275 for jj = 1:obj.n 272 for jj = 1:obj.n
276 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); 273 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
277 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); 274 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
278 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); 275 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
279 end 276 end
280 end 277 end
281 278
282 D = sparse(Dret); 279 D = sparse(Dret);
283 V = sparse(Vret); 280 V = sparse(Vret);
284 Vi = sparse(Viret); 281 Vi = sparse(Viret);
285 V = obj.evaluateCoefficientMatrix(V,x,y); 282 V = obj.evaluateCoefficientMatrix(V,x,y);
286 Vi = obj.evaluateCoefficientMatrix(Vi,x,y); 283 Vi = obj.evaluateCoefficientMatrix(Vi,x,y);
287 D = obj.evaluateCoefficientMatrix(D,x,y); 284 D = obj.evaluateCoefficientMatrix(D,x,y);
288 DD = diag(D); 285 DD = diag(D);
289 286
290 poseig = (DD>0); 287 poseig = (DD>0);
291 zeroeig = (DD==0); 288 zeroeig = (DD==0);
292 negeig = (DD<0); 289 negeig = (DD<0);
293 290
294 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); 291 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
295 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; 292 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
296 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 293 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
297 signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; 294 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
298 end 295 end
299 296
297 % Returns the boundary operator op for the boundary specified by the string boundary.
298 % op -- string or a cell array of strings
299 % boundary -- string
300 function varargout = getBoundaryOperator(obj, op, boundary)
301 assertIsMember(boundary, {'w', 'e', 's', 'n'})
302
303 if ~iscell(op)
304 op = {op};
305 end
306
307 for i = 1:numel(op)
308 switch op{i}
309 case 'e'
310 switch boundary
311 case 'w'
312 e = obj.e_w;
313 case 'e'
314 e = obj.e_e;
315 case 's'
316 e = obj.e_s;
317 case 'n'
318 e = obj.e_n;
319 end
320 varargout{i} = e;
321 end
322 end
323 end
324
325 % Returns square boundary quadrature matrix, of dimension
326 % corresponding to the number of boundary points
327 %
328 % boundary -- string
329 function H_b = getBoundaryQuadrature(obj, boundary)
330 assertIsMember(boundary, {'w', 'e', 's', 'n'})
331
332 e = obj.getBoundaryOperator('e', boundary);
333
334 switch boundary
335 case 'w'
336 H_b = inv(e'*obj.Hyi*e);
337 case 'e'
338 H_b = inv(e'*obj.Hyi*e);
339 case 's'
340 H_b = inv(e'*obj.Hxi*e);
341 case 'n'
342 H_b = inv(e'*obj.Hxi*e);
343 end
344 end
345
300 end 346 end
301 end 347 end