Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2d.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 Hypsyst2d < 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 % Grid
7 X,Y % Values of x and y for each grid point
8 order % Order accuracy for the approximation
9
10 D % non-stabalized scheme operator
11 A, B, E %Coefficient matrices
12
13 H % Discrete norm
14 % Norms in the x and y directions
15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
16 I_x,I_y, I_N
17 e_w, e_e, e_s, e_n
18 params %parameters for the coeficient matrice
19 end
20
21 methods
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)
24 default_arg('E', [])
25 xlim = lim{1};
26 ylim = lim{2};
27
28 if length(m) == 1
29 m = [m m];
30 end
31
32 obj.A=A;
33 obj.B=B;
34 obj.E=E;
35
36 m_x = m(1);
37 m_y = m(2);
38 obj.params = params;
39
40 ops_x = sbp.D2Standard(m_x,xlim,order);
41 ops_y = sbp.D2Standard(m_y,ylim,order);
42
43 obj.x = ops_x.x;
44 obj.y = ops_y.x;
45
46 obj.X = kr(obj.x,ones(m_y,1));
47 obj.Y = kr(ones(m_x,1),obj.y);
48
49 Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y);
50 Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y);
51 Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y);
52
53 obj.n = length(A(obj.params,0,0));
54
55 I_n = eye(obj.n);I_x = speye(m_x);
56 obj.I_x = I_x;
57 I_y = speye(m_y);
58 obj.I_y = I_y;
59
60
61 D1_x = kr(I_n, ops_x.D1, I_y);
62 obj.Hxi = kr(I_n, ops_x.HI, I_y);
63 D1_y = kr(I_n, I_x, ops_y.D1);
64 obj.Hyi = kr(I_n, I_x, ops_y.HI);
65
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);
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);
70
71 obj.m = m;
72 obj.h = [ops_x.h ops_y.h];
73 obj.order = order;
74
75 obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated;
76
77 end
78
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.
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.
83 % data is a function returning the data that should be applied at the boundary.
84 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
85 default_arg('type','char');
86 switch type
87 case{'c','char'}
88 [closure,penalty] = boundary_condition_char(obj,boundary);
89 case{'general'}
90 [closure,penalty] = boundary_condition_general(obj,boundary,L);
91 otherwise
92 error('No such boundary condition')
93 end
94 end
95
96 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
97 error('An interface function does not exist yet');
98 end
99
100 function N = size(obj)
101 N = obj.m;
102 end
103
104 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y)
105 params = obj.params;
106
107 if isa(mat,'function_handle')
108 [rows,cols] = size(mat(params,0,0));
109 matVec = mat(params,X',Y');
110 matVec = sparse(matVec);
111 side = max(length(X),length(Y));
112 else
113 matVec = mat;
114 [rows,cols] = size(matVec);
115 side = max(length(X),length(Y));
116 cols = cols/side;
117 end
118 ret = cell(rows,cols);
119
120 for ii = 1:rows
121 for jj=1:cols
122 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
123 end
124 end
125 ret = cell2mat(ret);
126 end
127
128 %Characteristic boundary conditions
129 function [closure, penalty] = boundary_condition_char(obj,boundary)
130 params = obj.params;
131 x = obj.x;
132 y = obj.y;
133
134 switch boundary
135 case {'w','W','west'}
136 e_ = obj.e_w;
137 mat = obj.A;
138 boundPos = 'l';
139 Hi = obj.Hxi;
140 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
141 side = max(length(y));
142 case {'e','E','east'}
143 e_ = obj.e_e;
144 mat = obj.A;
145 boundPos = 'r';
146 Hi = obj.Hxi;
147 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
148 side = max(length(y));
149 case {'s','S','south'}
150 e_ = obj.e_s;
151 mat = obj.B;
152 boundPos = 'l';
153 Hi = obj.Hyi;
154 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
155 side = max(length(x));
156 case {'n','N','north'}
157 e_ = obj.e_n;
158 mat = obj.B;
159 boundPos = 'r';
160 Hi = obj.Hyi;
161 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
162 side = max(length(x));
163 end
164 pos = signVec(1);
165 zeroval = signVec(2);
166 neg = signVec(3);
167
168 switch boundPos
169 case {'l'}
170 tau = sparse(obj.n*side,pos);
171 Vi_plus = Vi(1:pos,:);
172 tau(1:pos,:) = -abs(D(1:pos,1:pos));
173 closure = Hi*e_*V*tau*Vi_plus*e_';
174 penalty = -Hi*e_*V*tau*Vi_plus;
175 case {'r'}
176 tau = sparse(obj.n*side,neg);
177 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(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_';
180 penalty = -Hi*e_*V*tau*Vi_minus;
181 end
182 end
183
184 % General boundary condition in the form Lu=g(x)
185 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
186 params = obj.params;
187 x = obj.x;
188 y = obj.y;
189
190 switch boundary
191 case {'w','W','west'}
192 e_ = obj.e_w;
193 mat = obj.A;
194 boundPos = 'l';
195 Hi = obj.Hxi;
196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
197 L = obj.evaluateCoefficientMatrix(L,x(1),y);
198 side = max(length(y));
199 case {'e','E','east'}
200 e_ = obj.e_e;
201 mat = obj.A;
202 boundPos = 'r';
203 Hi = obj.Hxi;
204 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
205 L = obj.evaluateCoefficientMatrix(L,x(end),y);
206 side = max(length(y));
207 case {'s','S','south'}
208 e_ = obj.e_s;
209 mat = obj.B;
210 boundPos = 'l';
211 Hi = obj.Hyi;
212 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
213 L = obj.evaluateCoefficientMatrix(L,x,y(1));
214 side = max(length(x));
215 case {'n','N','north'}
216 e_ = obj.e_n;
217 mat = obj.B;
218 boundPos = 'r';
219 Hi = obj.Hyi;
220 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
221 L = obj.evaluateCoefficientMatrix(L,x,y(end));
222 side = max(length(x));
223 end
224
225 pos = signVec(1);
226 zeroval = signVec(2);
227 neg = signVec(3);
228
229 switch boundPos
230 case {'l'}
231 tau = sparse(obj.n*side,pos);
232 Vi_plus = Vi(1:pos,:);
233 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
234 V_plus = V(:,1:pos);
235 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
236
237 tau(1:pos,:) = -abs(D(1:pos,1:pos));
238 R = -inv(L*V_plus)*(L*V_minus);
239 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
240 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
241 case {'r'}
242 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));
244 Vi_plus = Vi(1:pos,:);
245 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
246
247 V_plus = V(:,1:pos);
248 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
249 R = -inv(L*V_minus)*(L*V_plus);
250 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
251 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
252 end
253 end
254
255 % 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
257 % [d+ ]
258 % D = [ d0 ]
259 % [ d-]
260 % 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)
262 params = obj.params;
263 syms xs ys
264 [V, D]= eig(mat(params,xs,ys));
265 Vi = inv(V);
266 xs = x;
267 ys = y;
268
269 side = max(length(x),length(y));
270 Dret = zeros(obj.n,side*obj.n);
271 Vret = zeros(obj.n,side*obj.n);
272 Viret = zeros(obj.n,side*obj.n);
273
274 for ii = 1:obj.n
275 for jj = 1:obj.n
276 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));
278 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
279 end
280 end
281
282 D = sparse(Dret);
283 V = sparse(Vret);
284 Vi = sparse(Viret);
285 V = obj.evaluateCoefficientMatrix(V,x,y);
286 Vi = obj.evaluateCoefficientMatrix(Vi,x,y);
287 D = obj.evaluateCoefficientMatrix(D,x,y);
288 DD = diag(D);
289
290 poseig = (DD>0);
291 zeroeig = (DD==0);
292 negeig = (DD<0);
293
294 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
295 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
296 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
297 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
298 end
299
300 end
301 end