Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2d.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 9b3d7fc61a36
children 459eeb99130f
comparison
equal deleted inserted replaced
368:53abf04f5e4e 369:9d1fc984f40d
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 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
22 methods 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) 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.
84 function [closure, penalty] = boundary_condition(obj,boundary,type,L) 84 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
85 default_arg('type','char'); 85 default_arg('type','char');
86 switch type 86 switch type
87 case{'c','char'} 87 case{'c','char'}
88 [closure,penalty]=boundary_condition_char(obj,boundary); 88 [closure,penalty] = boundary_condition_char(obj,boundary);
89 case{'general'} 89 case{'general'}
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)
97 error('An interface function does not exist yet'); 97 error('An interface function does not exist yet');
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));
112 else 112 else
113 matVec=mat; 113 matVec = mat;
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 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; y=obj.y; 131 x = obj.x;
132 132 y = obj.y;
133
133 switch boundary 134 switch boundary
134 case {'w','W','west'} 135 case {'w','W','west'}
135 e_=obj.e_w; 136 e_ = obj.e_w;
136 mat=obj.A; 137 mat = obj.A;
137 boundPos='l'; 138 boundPos = 'l';
138 Hi=obj.Hxi; 139 Hi = obj.Hxi;
139 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 140 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
140 side=max(length(y)); 141 side = max(length(y));
141 case {'e','E','east'} 142 case {'e','E','east'}
142 e_=obj.e_e; 143 e_ = obj.e_e;
143 mat=obj.A; 144 mat = obj.A;
144 boundPos='r'; 145 boundPos = 'r';
145 Hi=obj.Hxi; 146 Hi = obj.Hxi;
146 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 147 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
147 side=max(length(y)); 148 side = max(length(y));
148 case {'s','S','south'} 149 case {'s','S','south'}
149 e_=obj.e_s; 150 e_ = obj.e_s;
150 mat=obj.B; 151 mat = obj.B;
151 boundPos='l'; 152 boundPos = 'l';
152 Hi=obj.Hyi; 153 Hi = obj.Hyi;
153 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 154 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
154 side=max(length(x)); 155 side = max(length(x));
155 case {'n','N','north'} 156 case {'n','N','north'}
156 e_=obj.e_n; 157 e_ = obj.e_n;
157 mat=obj.B; 158 mat = obj.B;
158 boundPos='r'; 159 boundPos = 'r';
159 Hi=obj.Hyi; 160 Hi = obj.Hyi;
160 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 161 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
161 side=max(length(x)); 162 side = max(length(x));
162 end 163 end
163 164 pos = signVec(1);
164 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 165 zeroval = signVec(2);
165 166 neg = signVec(3);
167
166 switch boundPos 168 switch boundPos
167 case {'l'} 169 case {'l'}
168 tau=sparse(obj.n*side,pos); 170 tau = sparse(obj.n*side,pos);
169 Vi_plus=Vi(1:pos,:); 171 Vi_plus = Vi(1:pos,:);
170 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 172 tau(1:pos,:) = -abs(D(1:pos,1:pos));
171 closure=Hi*e_*V*tau*Vi_plus*e_'; 173 closure = Hi*e_*V*tau*Vi_plus*e_';
172 penalty=-Hi*e_*V*tau*Vi_plus; 174 penalty = -Hi*e_*V*tau*Vi_plus;
173 case {'r'} 175 case {'r'}
174 tau=sparse(obj.n*side,neg); 176 tau = sparse(obj.n*side,neg);
175 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 177 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
176 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 178 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
177 closure=Hi*e_*V*tau*Vi_minus*e_'; 179 closure = Hi*e_*V*tau*Vi_minus*e_';
178 penalty=-Hi*e_*V*tau*Vi_minus; 180 penalty = -Hi*e_*V*tau*Vi_minus;
179 end 181 end
180 end 182 end
181 183
182 184 % General boundary condition in the form Lu=g(x)
183 function [closure,penalty]=boundary_condition_general(obj,boundary,L) 185 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
184 params=obj.params; 186 params = obj.params;
185 x=obj.x; y=obj.y; 187 x = obj.x;
186 188 y = obj.y;
189
187 switch boundary 190 switch boundary
188 case {'w','W','west'} 191 case {'w','W','west'}
189 e_=obj.e_w; 192 e_ = obj.e_w;
190 mat=obj.A; 193 mat = obj.A;
191 boundPos='l'; 194 boundPos = 'l';
192 Hi=obj.Hxi; 195 Hi = obj.Hxi;
193 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
194 L=obj.evaluateCoefficientMatrix(L,x(1),y); 197 L = obj.evaluateCoefficientMatrix(L,x(1),y);
195 side=max(length(y)); 198 side = max(length(y));
196 case {'e','E','east'} 199 case {'e','E','east'}
197 e_=obj.e_e; 200 e_ = obj.e_e;
198 mat=obj.A; 201 mat = obj.A;
199 boundPos='r'; 202 boundPos = 'r';
200 Hi=obj.Hxi; 203 Hi = obj.Hxi;
201 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 204 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
202 L=obj.evaluateCoefficientMatrix(L,x(end),y); 205 L = obj.evaluateCoefficientMatrix(L,x(end),y);
203 side=max(length(y)); 206 side = max(length(y));
204 case {'s','S','south'} 207 case {'s','S','south'}
205 e_=obj.e_s; 208 e_ = obj.e_s;
206 mat=obj.B; 209 mat = obj.B;
207 boundPos='l'; 210 boundPos = 'l';
208 Hi=obj.Hyi; 211 Hi = obj.Hyi;
209 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 212 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
210 L=obj.evaluateCoefficientMatrix(L,x,y(1)); 213 L = obj.evaluateCoefficientMatrix(L,x,y(1));
211 side=max(length(x)); 214 side = max(length(x));
212 case {'n','N','north'} 215 case {'n','N','north'}
213 e_=obj.e_n; 216 e_ = obj.e_n;
214 mat=obj.B; 217 mat = obj.B;
215 boundPos='r'; 218 boundPos = 'r';
216 Hi=obj.Hyi; 219 Hi = obj.Hyi;
217 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 220 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
218 L=obj.evaluateCoefficientMatrix(L,x,y(end)); 221 L = obj.evaluateCoefficientMatrix(L,x,y(end));
219 side=max(length(x)); 222 side = max(length(x));
220 end 223 end
221 224
222 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 225 pos = signVec(1);
223 226 zeroval = signVec(2);
227 neg = signVec(3);
228
224 switch boundPos 229 switch boundPos
225 case {'l'} 230 case {'l'}
226 tau=sparse(obj.n*side,pos); 231 tau = sparse(obj.n*side,pos);
227 Vi_plus=Vi(1:pos,:); 232 Vi_plus = Vi(1:pos,:);
228 Vi_minus=Vi(pos+zeroval+1:obj.n*side,:); 233 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
229 V_plus=V(:,1:pos); 234 V_plus = V(:,1:pos);
230 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 235 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
231 236
232 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 237 tau(1:pos,:) = -abs(D(1:pos,1:pos));
233 R=-inv(L*V_plus)*(L*V_minus); 238 R = -inv(L*V_plus)*(L*V_minus);
234 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 239 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
235 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 240 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
236 case {'r'} 241 case {'r'}
237 tau=sparse(obj.n*side,neg); 242 tau = sparse(obj.n*side,neg);
238 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 243 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
239 Vi_plus=Vi(1:pos,:); 244 Vi_plus = Vi(1:pos,:);
240 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 245 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
241 246
242 V_plus=V(:,1:pos); 247 V_plus = V(:,1:pos);
243 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 248 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
244 R=-inv(L*V_minus)*(L*V_plus); 249 R = -inv(L*V_minus)*(L*V_plus);
245 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 250 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
246 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 251 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
247 end 252 end
248 end 253 end
249 254
250 255 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
251 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y) 256 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
252 params=obj.params; 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;
253 syms xs ys 263 syms xs ys
254 [V, D]=eig(mat(params,xs,ys)); 264 [V, D]= eig(mat(params,xs,ys));
255 Vi=inv(V); 265 Vi = inv(V);
256 xs=x; 266 xs = x;
257 ys=y; 267 ys = y;
258 268
259 side=max(length(x),length(y)); 269 side = max(length(x),length(y));
260 Dret=zeros(obj.n,side*obj.n); 270 Dret = zeros(obj.n,side*obj.n);
261 Vret=zeros(obj.n,side*obj.n); 271 Vret = zeros(obj.n,side*obj.n);
262 Viret=zeros(obj.n,side*obj.n); 272 Viret = zeros(obj.n,side*obj.n);
263 for ii=1:obj.n 273
264 for jj=1:obj.n 274 for ii = 1:obj.n
265 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 275 for jj = 1:obj.n
266 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 276 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
267 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(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));
268 end 279 end
269 end 280 end
270 281
271 D=sparse(Dret); 282 D = sparse(Dret);
272 V=sparse(Vret); 283 V = sparse(Vret);
273 Vi=sparse(Viret); 284 Vi = sparse(Viret);
274 V=obj.evaluateCoefficientMatrix(V,x,y); 285 V = obj.evaluateCoefficientMatrix(V,x,y);
275 Vi=obj.evaluateCoefficientMatrix(Vi,x,y); 286 Vi = obj.evaluateCoefficientMatrix(Vi,x,y);
276 D=obj.evaluateCoefficientMatrix(D,x,y); 287 D = obj.evaluateCoefficientMatrix(D,x,y);
277 DD=diag(D); 288 DD = diag(D);
278 289
279 poseig=(DD>0); 290 poseig = (DD>0);
280 zeroeig=(DD==0); 291 zeroeig = (DD==0);
281 negeig=(DD<0); 292 negeig = (DD<0);
282 293
283 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); 294 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
284 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 295 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
285 Vi=[Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 296 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
286 signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; 297 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
287 end 298 end
288 299
289 end 300 end
290 end 301 end