Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 59:e431c1260f52

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Euler1d: Changed penalties to be calcualted numerically.
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 14 Nov 2015 11:59:49 -0800
parents 24103284e09d
children e707cd9e6d0a
comparison
equal deleted inserted replaced
58:24103284e09d 59:e431c1260f52
232 pt(p) = 1:length(p); % Inverse permutation 232 pt(p) = 1:length(p); % Inverse permutation
233 233
234 gamma = obj.gamma; 234 gamma = obj.gamma;
235 235
236 function o = closure_fun(q,t) 236 function o = closure_fun(q,t)
237 % Extract solution at the boundary
237 q_s = e_S'*q; 238 q_s = e_S'*q;
238 rho = q_s(1); 239 rho = q_s(1);
239 u = q_s(2)/rho; 240 u = q_s(2)/rho;
240 e = q_s(3); 241 e = q_s(3);
241 242
242 c = obj.c(q_s); 243 c = obj.c(q_s);
243 244
245 % Calculate transformation matrix
244 T = obj.T(q_s); 246 T = obj.T(q_s);
245 247 Tin = T(:,p_in);
246 % L = [ 248 Tot = T(:,p_ot);
247 % 0 1 0; 249
248 % 0 -1/2*u 1; 250 % Convert bc from ordinary form to characteristic form.
249 % ]; 251 % Lq = g => w_in = Rw_ot + g_tilde
250 LTpi = [ 252 L = [
251 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)); 253 0 1 0;
252 0, 2/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2)); 254 0 -1/2*u 1;
253 ]; 255 ];
254 R = [ 256 g = [
255 (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)); 257 v_data(t);
256 (2*c*rho*s*(u - 2))/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2)) - 1; 258 p_data(t);
257 ]; 259 ];
258 260
261 R =-inv(L*Tin)*(L*Tot);
262 g_tilde = inv(L*Tin)*g;
263
264 % Setup the penalty parameter
259 tau1 = -2*[ 265 tau1 = -2*[
260 u, 0; 266 u, 0;
261 0, c + abs(u); 267 0, c + abs(u);
262 ]; 268 ];
263 tau2 = [0 0]; % Penalty only on ingoing char. 269 tau2 = [0 0]; % Penalty only on ingoing char.
264 270
265
266 tauHat = [tau1; tau2]; 271 tauHat = [tau1; tau2];
267
268 tau = -s*e_S*sparse(T*tauHat(pt,:)); 272 tau = -s*e_S*sparse(T*tauHat(pt,:));
269 273
274 % Calculate the numerical caracteristics and apply the penalty
270 w_s = inv(T)*q_s; 275 w_s = inv(T)*q_s;
271 % w_s = T\q_s; 276 % w_s = T\q_s;
272 % w_s = Tinv * q_s; % Med analytisk matris 277 % w_s = Tinv * q_s; % Med analytisk matris
273 w_in = w_s(p_in); 278 w_in = w_s(p_in);
274 w_ot = w_s(p_ot); 279 w_ot = w_s(p_ot);
275 280
276 281 o = obj.Hi * tau * (w_in - R*w_ot - g_tilde);
277 o = obj.Hi * tau * (w_in - R*w_ot - LTpi*[v_data(t); p_data(t)]);
278 end 282 end
279 283
280 closure = @closure_fun; 284 closure = @closure_fun;
281 end 285 end
282 286
298 p = [p_in, p_ot]; % Permutation to sort 302 p = [p_in, p_ot]; % Permutation to sort
299 pt(p) = 1:length(p); % Inverse permutation 303 pt(p) = 1:length(p); % Inverse permutation
300 304
301 gamma = obj.gamma; 305 gamma = obj.gamma;
302 306
307
303 function o = closure_fun(q,t) 308 function o = closure_fun(q,t)
309 % Extract solution at the boundary
304 q_s = e_S'*q; 310 q_s = e_S'*q;
305 rho = q_s(1); 311 rho = q_s(1);
306 u = q_s(2)/rho; 312 u = q_s(2)/rho;
307 e = q_s(3); 313 e = q_s(3);
308 314
309 c = obj.c(q_s); 315 c = obj.c(q_s);
310 316
311 317 % Calculate transformation matrix
312 T = obj.T(q_s); 318 T = obj.T(q_s);
313 319 Tin = T(:,p_in);
320 Tot = T(:,p_ot);
321
322 % Convert bc from ordinary form to characteristic form.
323 % Lq = g => w_in = Rw_ot + g_tilde
314 L = (gamma -1)*[0 -1/2*u 1]; 324 L = (gamma -1)*[0 -1/2*u 1];
315 LTpi = 1/(L*T(:,p_in)); 325 g = [p_data(t)];
316 R = -[0 (gamma*(e - rho*u^2/2) - rho*u*c)/(gamma*(e - rho*u^2/2) + rho*u*c)]; 326
317 327 R =-inv(L*Tin)*(L*Tot);
328 g_tilde = inv(L*Tin)*g;
329
330 % Setup the penalty parameter
318 tau1 = [-2*(c - abs(u))]; 331 tau1 = [-2*(c - abs(u))];
319 tau2 = [0; 0]; % Penalty only on ingoing char. 332 tau2 = [0; 0]; % Penalty only on ingoing char.
320 333
321
322 tauHat = [tau1; tau2]; 334 tauHat = [tau1; tau2];
323 tau = -s*e_S*sparse(T*tauHat(pt)); 335 tau = -s*e_S*sparse(T*tauHat(pt,:));
324 336
337 % Calculate the numerical caracteristics and apply the penalty
325 w_s = inv(T)*q_s; 338 w_s = inv(T)*q_s;
326 % w_s = T\q_s; 339 % w_s = T\q_s;
327 % w_s = Tinv * q_s; % Med analytisk matris 340 % w_s = Tinv * q_s; % Med analytisk matris
328 w_in = w_s(p_in); 341 w_in = w_s(p_in);
329 w_ot = w_s(p_ot); 342 w_ot = w_s(p_ot);
330 343
331 o = obj.Hi * tau * (w_in - R*w_ot - LTpi*p_data(t)); 344 o = obj.Hi * tau * (w_in - R*w_ot - g_tilde);
332 end 345 end
333 346
334 closure = @closure_fun; 347 closure = @closure_fun;
335 348
336 end 349 end