Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 54:2194cd385419

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 11 Nov 2015 19:21:38 -0800
parents 75ebf5d3cfe5
children 9a647dcccbdd
comparison
equal deleted inserted replaced
53:cca09cc5121d 54:2194cd385419
31 gamma = gama; 31 gamma = gama;
32 32
33 [x, h] = util.get_grid(xlim{:},m); 33 [x, h] = util.get_grid(xlim{:},m);
34 34
35 if do_upwind 35 if do_upwind
36 ops = spb.Upwind(m,h,order); 36 ops = sbp.Upwind(m,h,order);
37 Dp = ops.derivatives.Dp; 37 Dp = ops.derivatives.Dp;
38 Dm = ops.derivatives.Dm; 38 Dm = ops.derivatives.Dm;
39 39
40 printExpr('issparse(Dp)');
41 printExpr('issparse(Dm)');
42
43 D1 = (Dp + Dm)/2; 40 D1 = (Dp + Dm)/2;
41 Ddisp = (Dp - Dm)/2;
44 else 42 else
45 ops = opsGen(m,h,order); 43 ops = opsGen(m,h,order);
46 D1 = sparse(ops.derivatives.D1); 44 D1 = sparse(ops.derivatives.D1);
47 end 45 end
48 46
51 e_l = sparse(ops.boundary.e_1); 49 e_l = sparse(ops.boundary.e_1);
52 e_r = sparse(ops.boundary.e_m); 50 e_r = sparse(ops.boundary.e_m);
53 51
54 I_x = speye(m); 52 I_x = speye(m);
55 I_3 = speye(3); 53 I_3 = speye(3);
56
57 54
58 D1 = kr(D1, I_3); 55 D1 = kr(D1, I_3);
59 if do_upwind 56 if do_upwind
60 Ddisp = kr(Ddisp,I_3); 57 Ddisp = kr(Ddisp,I_3);
61 end 58 end
94 function o = c(Q) 91 function o = c(Q)
95 % Speed of light c = sqrt(obj.gamma*p/rho); 92 % Speed of light c = sqrt(obj.gamma*p/rho);
96 o = sqrt(gamma*p(Q)./Q(1,:)); 93 o = sqrt(gamma*p(Q)./Q(1,:));
97 end 94 end
98 95
99
100 % R =
101 % [sqrt(2*(gamma-1))*rho , rho , rho ;
102 % sqrt(2*(gamma-1))*rho*u , rho*(u+c) , rho*(u-c) ;
103 % sqrt(2*(gamma-1))*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c, e+(gamma-1)*(e-rho*u^2/2)-rho*u*c]);
104 function o = T(q) 96 function o = T(q)
97 % T is the transformation such that A = T*Lambda*inv(T)
98 % where Lambda=diag(u, u+c, u-c)
105 rho = q(1); 99 rho = q(1);
106 u = q(2)/q(1); 100 u = q(2)/q(1);
107 e = q(3); 101 e = q(3);
108 102
109 c = sqrt(gamma*p(q)/rho); 103 c = sqrt(gamma*p(q)/rho);
111 sqrt2gamm = sqrt(2*(gamma-1)); 105 sqrt2gamm = sqrt(2*(gamma-1));
112 106
113 o = [ 107 o = [
114 sqrt2gamm*rho , rho , rho ; 108 sqrt2gamm*rho , rho , rho ;
115 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ; 109 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ;
116 sqrt2gamm*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c , e+(gamma-1)*(e-rho*u^2/2)-rho*u*c 110 sqrt2gamm*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c , e+(gamma-1)*(e-rho*u^2/2)-rho*u*c ;
117 ]; 111 ];
118 % Devide columns by stuff to get rid of extra comp? 112 % Devide columns by stuff to get rid of extra comp?
119 end 113 end
120 114
121 function o = Fx(q) 115 function o = Fx(q)
122 Q = reshape(q,3,m); 116 Q = reshape(q,3,m);
123 o = reshape(F(Q),3*m,1); 117 o = reshape(F(Q),3*m,1);
124 o = D1*o; 118 o = D1*o;
125 end 119 end
126 120
127 121 function o = Fx_disp(q)
128 % A=R*Lambda*inv(R), där Lambda=diag(u, u+c, u-c) (c är ljudhastigheten) 122 Q = reshape(q,3,m);
129 % c^2=gamma*p/rho 123 f = reshape(F(Q),3*m,1);
130 % function o = A(rho,u,e) 124
131 % end 125 c = c(Q);
126 lambda_max = c+abs(Q(2,:)./Q(1,:));
127 % lambda_max = max(lambda_max);
128
129 lamb_Q(1,:) = lambda_max.*Q(1,:);
130 lamb_Q(2,:) = lambda_max.*Q(2,:);
131 lamb_Q(3,:) = lambda_max.*Q(3,:);
132
133 lamb_q = reshape(lamb_Q,3*m, 1);
134
135 o = -D1*f + Ddisp*lamb_q;
136 end
137
132 if do_upwind 138 if do_upwind
133 obj.D = @(q)-Fx(q) + Ddisp*(1)*u; 139 obj.D = @Fx_disp;
134 else 140 else
135 obj.D = @(q)-Fx(q); 141 obj.D = @(q)-Fx(q);
136 end 142 end
137 143
138 obj.Fx = @Fx; 144 obj.Fx = @Fx;
151 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 157 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
152 158
153 % Boundary condition on form 159 % Boundary condition on form
154 % w_in = R*w_out + g, where g is data 160 % w_in = R*w_out + g, where g is data
155 161
156 % How to handle when the number of BC we want to set changes
157 % How to handel u = 0 as an example
158
159 % Måste sätta in s och fixa tecken som de ska vara
160
161 % Kanske ska man tala om vilka karakteristikor som är in och ut i anropet
162 % Man kan sen kolla om det stämmer (men hur blir det med u=0?)
163 % Och antingen tillåta att man skickar in flera alternativ och väljer automatiskt
164 % eller låta koden utanför bestämma vilken penalty som ska appliceras.
165
166 switch type 162 switch type
167 case 'wall' 163 case 'wall'
168 closure = obj.boundary_condition_wall(boundary); 164 closure = obj.boundary_condition_wall(boundary);
165 case 'inflow'
166 closure = obj.boundary_condition_inflow(boundary,varargin{:});
167 case 'outflow'
168 closure = obj.boundary_condition_outflow(boundary,varargin{:});
169 otherwise 169 otherwise
170 error('Unsupported bc type: %s', type); 170 error('Unsupported bc type: %s', type);
171 end 171 end
172 172
173 % T = obj.T(e_S*v); 173 end
174 174
175 % c = sqrt(gamma*p/rho); 175 function closure = boundary_condition_inflow(obj, boundary, p_data, v_data)
176 % l = [u, u+c, u-c]; 176 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
177 177
178 % p_in = find(s*l <= 0); 178 switch s
179 % p_ot = find(s*l > 0); 179 case -1
180 180 p_in = [1 2];
181 % p = [p_in, p_ot] % Permutation to sort 181 p_ot = [3];
182 % pt(p) = 1:length(p); % Inverse permutation 182 case 1
183 183 p_in = [1 3];
184 % L_in = diag(abs(l(p_in))); 184 p_ot = [2];
185 % L_ot = diag(abs(l(p_ot))); 185 otherwise
186 186 error();
187 % Lambda = diag(u, u+c, u-c); 187 end
188 188
189 % tau = -e_S*sparse(T*[L_in; R'*L_in]); % Något med pt 189 p = [p_in, p_ot]; % Permutation to sort
190 190 pt(p) = 1:length(p); % Inverse permutation
191 191
192 % w_s = T'*(e_S*v); 192 gamma = obj.gamma;
193 % w_in = w_s(p_in); 193
194 % w_ot = w_s(p_ot); 194 function o = closure_fun(q,t)
195 195 q_s = e_S'*q;
196 196 rho = q_s(1);
197 % function closure_fun(q,t) 197 u = q_s(2)/rho;
198 % obj.Hi * tau * (w_in - R*w_ot - g(t)); 198 e = q_s(3);
199 % end 199
200 200 c = obj.c(q_s);
201 % function closure_fun_indep(q,t) 201
202 % obj.Hi * tau * (w_in - R*w_ot - g; 202 T = obj.T(q_s);
203 % end 203
204 204 % L = [
205 205 % 0 1 0;
206 % switch class(g) 206 % 0 -1/2*u 1;
207 % case 'double' 207 % ];
208 % closure = @closure_fun; 208 LTpi = [
209 % case 'function_handle' 209 2^(1/2)/(2*rho*u*(gamma - 1)^(1/2)), -(2^(1/2)*(u - c*s))/(u*(2*e*gamma - rho*(gamma*u^2 - c*s*u + 2*c*s))*(gamma - 1)^(1/2));
210 % closure = @closure_fun_indep; 210 0, 2/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2));
211 % otherwise 211 ];
212 % error('Wierd data argument!'); 212 R = [
213 % end 213 (2^(1/2)*c*s*(rho*(gamma - 1)*u^2 + 2*rho*u - 2*e*gamma))/(u*(gamma - 1)^(1/2)*(- gamma*rho*u^2 + c*rho*s*u + 2*e*gamma - 2*c*rho*s));
214 (2*c*rho*s*(u - 2))/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2)) - 1;
215 ];
216
217 tau1 = -2*[
218 u, 0;
219 0, c + abs(u);
220 ];
221 tau2 = [0 0]; % Penalty only on ingoing char.
222
223
224 tauHat = [tau1; tau2];
225 tau = -s*e_S*sparse(T*tauHat(pt,:));
226
227 w_s = inv(T)*q_s;
228 % w_s = T\q_s;
229 % w_s = Tinv * q_s; % Med analytisk matris
230 w_in = w_s(p_in);
231 w_ot = w_s(p_ot);
232
233
234 o = obj.Hi * tau * (w_in - R*w_ot - LTpi*[v_data(t); p_data(t)]);
235 end
236
237 closure = @closure_fun;
238 end
239
240
241 function closure = boundary_condition_outflow(obj, boundary, p_data)
242 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
243
244 switch s
245 case -1
246 p_in = 2;
247 p_ot = [1 3];
248 case 1
249 p_in = 3;
250 p_ot = [1 2];
251 otherwise
252 error();
253 end
254
255 p = [p_in, p_ot]; % Permutation to sort
256 pt(p) = 1:length(p); % Inverse permutation
257
258 gamma = obj.gamma;
259
260 function o = closure_fun(q,t)
261 q_s = e_S'*q;
262 rho = q_s(1);
263 u = q_s(2)/rho;
264 e = q_s(3);
265
266 c = obj.c(q_s);
267
268
269 T = obj.T(q_s);
270
271 L = (gamma -1)*[0 -1/2*u 1];
272 LTpi = 1/(L*T(:,p_in));
273 R = -[0 (gamma*(e - rho*u^2/2) - rho*u*c)/(gamma*(e - rho*u^2/2) + rho*u*c)];
274
275 tau1 = [-2*(c - abs(u))];
276 tau2 = [0; 0]; % Penalty only on ingoing char.
277
278
279 tauHat = [tau1; tau2];
280 tau = -s*e_S*sparse(T*tauHat(pt));
281
282 w_s = inv(T)*q_s;
283 % w_s = T\q_s;
284 % w_s = Tinv * q_s; % Med analytisk matris
285 w_in = w_s(p_in);
286 w_ot = w_s(p_ot);
287
288 o = obj.Hi * tau * (w_in - R*w_ot - LTpi*p_data(t));
289 end
290
291 closure = @closure_fun;
214 292
215 end 293 end
216 294
217 % Set wall boundary condition v = 0. 295 % Set wall boundary condition v = 0.
218 function closure = boundary_condition_wall(obj,boundary) 296 function closure = boundary_condition_wall(obj,boundary)
277 355
278 356
279 tau1 = -2*c; 357 tau1 = -2*c;
280 tau2 = [0; 0]; % Penalty only on ingoing char. 358 tau2 = [0; 0]; % Penalty only on ingoing char.
281 359
282 % L_in = diag(abs(l(p_in))); 360 % Lambda_in = diag(abs(l(p_in)));
283 % L_ot = diag(abs(l(p_ot))); 361 % Lambda_ot = diag(abs(l(p_ot)));
284 362
285 tauHat = [tau1; tau2]; 363 tauHat = [tau1; tau2];
286 tau = -s*e_S*sparse(T*tauHat(pt)); 364 tau = -s*e_S*sparse(T*tauHat(pt));
287 365
288 w_s = inv(T)*q_s; 366 w_s = inv(T)*q_s;