Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 423:a2cb0d4f4a02 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 07 Feb 2017 15:47:51 +0100
parents 9d1fc984f40d
children 459eeb99130f
comparison
equal deleted inserted replaced
218:da058ce66876 423:a2cb0d4f4a02
1 classdef Hypsyst2dCurve < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 n % size of system
5 h % Grid spacing
6 X,Y % Values of x and y for each grid point
7
8 J, Ji % Jacobaian and inverse Jacobian
9 xi,eta
10 Xi,Eta
11
12 A,B
13 X_eta, Y_eta
14 X_xi,Y_xi
15 order % Order accuracy for the approximation
16
17 D % non-stabalized scheme operator
18 Ahat, Bhat, E
19
20 H % Discrete norm
21 Hxii,Hetai % Kroneckerd norms in xi and eta.
22 I_xi,I_eta, I_N, onesN
23 e_w, e_e, e_s, e_n
24 index_w, index_e,index_s,index_n
25 params % Parameters for the coeficient matrice
26 end
27
28
29 methods
30 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
31 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti)
32 default_arg('E', [])
33 xilim = {0 1};
34 etalim = {0 1};
35
36 if length(m) == 1
37 m = [m m];
38 end
39 obj.params = params;
40 obj.A=A;
41 obj.B=B;
42
43 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta);
44 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi);
45 obj.E=@(params,x,y,~,~)E(params,x,y);
46
47 m_xi = m(1);
48 m_eta = m(2);
49 m_tot=m_xi*m_eta;
50
51 ops_xi = sbp.D2Standard(m_xi,xilim,order);
52 ops_eta = sbp.D2Standard(m_eta,etalim,order);
53
54 obj.xi = ops_xi.x;
55 obj.eta = ops_eta.x;
56
57 obj.Xi = kr(obj.xi,ones(m_eta,1));
58 obj.Eta = kr(ones(m_xi,1),obj.eta);
59
60 obj.n = length(A(obj.params,0,0));
61 obj.onesN=ones(obj.n);
62
63 obj.index_w=1:m_eta;
64 obj.index_e=(m_tot-m_e
65
66 metric_termsta+1):m_tot;
67 obj.index_s=1:m_eta:(m_tot-m_eta+1);
68 obj.index_n=(m_eta):m_eta:m_tot;
69
70 I_n = eye(obj.n);
71 I_xi = speye(m_xi);
72 obj.I_xi = I_xi;
73 I_eta = speye(m_eta);
74 obj.I_eta = I_eta;
75
76 D1_xi = kr(I_n, ops_xi.D1, I_eta);
77 obj.Hxii = kr(I_n, ops_xi.HI, I_eta);
78 D1_eta = kr(I_n, I_xi, ops_eta.D1);
79 obj.Hetai = kr(I_n, I_xi, ops_eta.HI);
80
81 obj.e_w = kr(I_n, ops_xi.e_l, I_eta);
82 obj.e_e = kr(I_n, ops_xi.e_r, I_eta);
83 obj.e_s = kr(I_n, I_xi, ops_eta.e_l);
84 obj.e_n = kr(I_n, I_xi,
85
86 metric_termsops_eta.e_r);
87
88 [X,Y] = ti.map(obj.xi,obj.eta);
89
90 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1);
91 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1);
92
93 obj.X = reshape(X,m_tot,1);
94 obj.Y = reshape(Y,m_tot,1);
95 obj.X_xi = reshape(x_xi,m_tot,1);
96 obj.Y_xi = reshape(y_xi,m_tot,1);
97 obj.X_eta = reshape(x_eta,m_tot,1);
98 obj.Y_eta = reshape(y_eta,m_tot,1);
99
100 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta);
101 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi);
102 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
103
104 obj.m = m;
105 obj.h = [ops_xi.h ops_eta.h];
106 obj.order = order;
107 obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi;
108 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
109
110 obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated;
111
112 end
113
114 % Closure functions return the opertors applied to the own doamin to close the boundary
115 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
116 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'.
117 % type is a string specifying the type of boundary condition if there are several.
118 % data is a function returning the data that should be applied at the boundary.
119 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
120 default_arg('type','char');
121 switch type
122 case{'c','char'}
123 [closure,penalty] = boundary_condition_char(obj,boundary);
124 case{'general'}
125 [closure,penalty] = boundary_condition_general(obj,boundary,L);
126 otherwise
127 error('No such boundary condition')
128 end
129 end
130
131 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions)
132 error('An interface function does not exist yet');
133 end
134
135 function N = size(obj)
136 N = obj.m;
137 end
138
139 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
140 params = obj.params;
141
142 if isa(mat,'function_handle')
143 [rows,cols] = size(mat(params,0,0,0,0));
144 x_ = kr(obj.onesN,x_);
145 y_ = kr(obj.onesN,y_);
146 matVec = mat(params,X',Y',x_',y_');
147 matVec = sparse(matVec);
148 side = max(length(X),length(Y));
149 else
150 matVec = mat;
151 [rows,cols] = size(matVec);
152 side = max(length(X),length(Y));
153 cols = cols/side;
154 end
155
156 ret = cell(rows,cols);
157 for ii = 1:rows
158 for jj = 1:cols
159 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
160 end
161 end
162 ret = cell2mat(ret);
163 end
164
165 %Characteristic boundary conditions
166 function [closure, penalty] = boundary_condition_char(obj,boundary)
167 params = obj.params;
168 X = obj.X;
169 Y = obj.Y;
170 xi = obj.xi;
171 eta = obj.eta;
172
173 switch boundary
174 case {'w','W','west'}
175 e_ = obj.e_w;
176 mat = obj.Ahat;
177 boundPos = 'l';
178 Hi = obj.Hxii;
179 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
180 side = max(length(eta));
181 case {'e','E','east'}
182 e_ = obj.e_e;
183 mat = obj.Ahat;
184 boundPos = 'r';
185 Hi = obj.Hxii;
186 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
187 side = max(length(eta));
188 case {'s','S','south'}
189 e_ = obj.e_s;
190 mat = obj.Bhat;
191 boundPos = 'l';
192 Hi = obj.Hetai;
193 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
194 side = max(length(xi));
195 case {'n','N','north'}
196 e_ = obj.e_n;
197 mat = obj.Bhat;
198 boundPos = 'r';
199 Hi = obj.Hetai;
200 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
201 side = max(length(xi));
202 end
203
204 pos = signVec(1);
205 zeroval = signVec(2);
206 neg = signVec(3);
207
208 switch boundPos
209 case {'l'}
210 tau = sparse(obj.n*side,pos);
211 Vi_plus = Vi(1:pos,:);
212 tau(1:pos,:) = -abs(D(1:pos,1:pos));
213 closure = Hi*e_*V*tau*Vi_plus*e_';
214 penalty = -Hi*e_*V*tau*Vi_plus;
215 case {'r'}
216 tau = sparse(obj.n*side,neg);
217 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
218 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
219 closure = Hi*e_*V*tau*Vi_minus*e_';
220 penalty = -Hi*e_*V*tau*Vi_minus;
221 end
222 end
223
224
225 % General boundary condition in the form Lu=g(x)
226 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
227 params = obj.params;
228 X = obj.X;
229 Y = obj.Y;
230 xi = obj.xi;
231 eta = obj.eta;
232
233 switch boundary
234 case {'w','W','west'}
235 e_ = obj.e_w;
236 mat = obj.Ahat;
237 boundPos = 'l';
238 Hi = obj.Hxii;
239 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
240
241 Ji_vec = diag(obj.Ji);
242 Ji = diag(Ji_vec(obj.index_w));
243 xi_x = Ji*obj.Y_eta(obj.index_w);
244 xi_y = -Ji*obj.X_eta(obj.index_w);
245 L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]);
246 side = max(length(eta));
247 case {'e','E','east'}
248 e_ = obj.e_e;
249 mat = obj.Ahat;
250 boundPos = 'r';
251 Hi = obj.Hxii;
252 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
253
254 Ji_vec = diag(obj.Ji);
255 Ji = diag(Ji_vec(obj.index_e));
256 xi_x = Ji*obj.Y_eta(obj.index_e);
257 xi_y = -Ji*obj.X_eta(obj.index_e);
258 L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]);
259 side = max(length(eta));
260 case {'s','S','south'}
261 e_ = obj.e_s;
262 mat = obj.Bhat;
263 boundPos = 'l';
264 Hi = obj.Hetai;
265 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
266
267 Ji_vec = diag(obj.Ji);
268 Ji = diag(Ji_vec(obj.index_s));
269 eta_x = Ji*obj.Y_xi(obj.index_s);
270 eta_y = -Ji*obj.X_xi(obj.index_s);
271 L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
272 side = max(length(xi));
273 case {'n','N','north'}
274 e_ = obj.e_n;
275 mat = obj.Bhat;
276 boundPos = 'r';
277 Hi = obj.Hetai;
278 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
279
280 Ji_vec = diag(obj.Ji);
281 Ji = diag(Ji_vec(obj.index_n));
282 eta_x = Ji*obj.Y_xi(obj.index_n);
283 eta_y = -Ji*obj.X_xi(obj.index_n);
284 L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]);
285 side = max(length(xi));
286 end
287
288 pos = signVec(1);
289 zeroval = signVec(2);
290 neg = signVec(3);
291
292 switch boundPos
293 case {'l'}
294 tau = sparse(obj.n*side,pos);
295 Vi_plus = Vi(1:pos,:);
296 Vi_minus = Vi(pos+1:obj.n*side,:);
297 V_plus = V(:,1:pos);
298 V_minus = V(:,(pos)+1:obj.n*side);
299
300 tau(1:pos,:) = -abs(D(1:pos,1:pos));
301 R = -inv(L*V_plus)*(L*V_minus);
302 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
303 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
304 case {'r'}
305 tau = sparse(obj.n*side,neg);
306 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
307 Vi_plus = Vi(1:pos,:);
308 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
309
310 V_plus = V(:,1:pos);
311 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
312 R = -inv(L*V_minus)*(L*V_plus);
313 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
314 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
315 end
316 end
317
318 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
319 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
320 % [d+ ]
321 % D = [ d0 ]
322 % [ d-]
323 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D
324 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y,x_,y_)
325 params = obj.params;
326 syms xs ys
327 if(sum(abs(x_)) ~= 0)
328 syms xs_
329 else
330 xs_ = 0;
331 end
332
333 if(sum(abs(y_))~= 0)
334 syms ys_;
335 else
336 ys_ = 0;
337 end
338
339 [V, D] = eig(mat(params,xs,ys,xs_,ys_));
340 Vi = inv(V);
341 syms xs ys xs_ ys_
342
343 xs = x;
344 ys = y;
345 xs_ = x_;
346 ys_ = y_;
347
348 side = max(length(x),length(y));
349 Dret = zeros(obj.n,side*obj.n);
350 Vret = zeros(obj.n,side*obj.n);
351 Viret = zeros(obj.n,side*obj.n);
352 for ii = 1:obj.n
353 for jj = 1:obj.n
354 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));
356 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
357 end
358 end
359
360 D = sparse(Dret);
361 V = sparse(Vret);
362 Vi = sparse(Viret);
363 V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
364 D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
365 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
366 DD = diag(D);
367
368 poseig = (DD>0);
369 zeroeig = (DD==0);
370 negeig = (DD<0);
371
372 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
373 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
374 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
375 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
376 end
377 end
378 end