Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 110:f5ed7ff58115

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 08 Dec 2015 13:03:04 +0100
parents ce4eecbcb915
children 0e66299592cc
comparison
equal deleted inserted replaced
109:142974097efc 110:f5ed7ff58115
224 224
225 end 225 end
226 226
227 227
228 % Sets the boundary condition Lq = g, where 228 % Sets the boundary condition Lq = g, where
229 % L = L(rho, u, e), g = g(t) 229 % L = L(rho, u, e)
230 % p_in are the indecies of the ingoing characteristics. 230 % p_in are the indecies of the ingoing characteristics.
231 function closure = boundary_condition_L(obj, boundary, L_fun, g_fun, p_in) 231 %
232 % Returns closure(q,g)
233 function closure = boundary_condition_L(obj, boundary, L_fun, p_in)
232 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 234 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
233 235
234 p_ot = 1:3; 236 p_ot = 1:3;
235 p_ot(p_in) = []; 237 p_ot(p_in) = [];
236 238
237 p = [p_in, p_ot]; % Permutation to sort 239 p = [p_in, p_ot]; % Permutation to sort
238 pt(p) = 1:length(p); % Inverse permutation 240 pt(p) = 1:length(p); % Inverse permutation
239 241
240 function o = closure_fun(q,t) 242 function o = closure_fun(q,g)
241 % Extract solution at the boundary 243 % Extract solution at the boundary
242 q_s = e_S'*q; 244 q_s = e_S'*q;
243 rho = q_s(1); 245 rho = q_s(1);
244 u = q_s(2)/rho; 246 u = q_s(2)/rho;
245 e = q_s(3); 247 e = q_s(3);
260 262
261 tauHat = [tau1; tau2]; 263 tauHat = [tau1; tau2];
262 tau = e_S*sparse(T*tauHat(pt,:)); 264 tau = e_S*sparse(T*tauHat(pt,:));
263 265
264 L = L_fun(rho,u,e); 266 L = L_fun(rho,u,e);
265 g = g_fun(t);
266 267
267 o = 1/2*obj.Hi * tau * inv(L*Tin)*(L*q_s - g); 268 o = 1/2*obj.Hi * tau * inv(L*Tin)*(L*q_s - g);
268 end 269 end
269 closure = @closure_fun; 270 closure = @closure_fun;
270 end 271 end
271 272
272 function closure = boundary_condition_char(obj,boundary,w_data) 273 % Return closure(q,g)
274 function closure = boundary_condition_char(obj,boundary)
273 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 275 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
274 276
275 function o = closure_fun(q,t) 277 function o = closure_fun(q, w_data)
276 q_s = e_S'*q; 278 q_s = e_S'*q;
277 rho = q_s(1); 279 rho = q_s(1);
278 u = q_s(2)/rho; 280 u = q_s(2)/rho;
279 e = q_s(3); 281 e = q_s(3);
280 282
298 tau = -s*e_S*sparse(T*tauHat(pt,:)); 300 tau = -s*e_S*sparse(T*tauHat(pt,:));
299 301
300 w_s = inv(T)*q_s; 302 w_s = inv(T)*q_s;
301 w_in = w_s(p_in); 303 w_in = w_s(p_in);
302 304
303 w_s_data = w_data(t); 305 w_in_data = w_data(p_in);
304 w_in_data = w_s_data(p_in);
305 306
306 o = 1/2*obj.Hi * tau * (w_in - w_in_data); 307 o = 1/2*obj.Hi * tau * (w_in - w_in_data);
307 end 308 end
308 309
309 closure = @closure_fun; 310 closure = @closure_fun;
310 311 end
311 end 312
312 313
313 function closure = boundary_condition_inflow(obj, boundary, p_data, v_data) 314 % Return closure(q,[v; p])
315 function closure = boundary_condition_inflow(obj, boundary)
314 [~,~,s] = obj.get_boundary_ops(boundary); 316 [~,~,s] = obj.get_boundary_ops(boundary);
315 317
316 switch s 318 switch s
317 case -1 319 case -1
318 p_in = [1 2]; 320 p_in = [1 2];
320 p_in = [1 3]; 322 p_in = [1 3];
321 end 323 end
322 324
323 a = obj.gamma - 1; 325 a = obj.gamma - 1;
324 L = @(rho,u,~)[ 326 L = @(rho,u,~)[
325 0 1/rho 0; 327 0 1/rho 0; %v
326 0 -1/2*u*a a; 328 0 -1/2*u*a a; %p
327 ]; 329 ];
328 g = @(t)[ 330
329 v_data(t); 331 closure_raw = boundary_condition_L(obj, boundary, L, g, p_in);
330 p_data(t); 332 closure = @(q,p,v) closure_raw(q,[v; p]);
331 ]; 333 end
332 334
333 closure = boundary_condition_L(obj, boundary, L, g, p_in); 335 % Return closure(q, p)
334 end 336 function closure = boundary_condition_outflow(obj, boundary)
335
336
337 function closure = boundary_condition_outflow(obj, boundary, p_data)
338 [~,~,s] = obj.get_boundary_ops(boundary); 337 [~,~,s] = obj.get_boundary_ops(boundary);
339 338
340 switch s 339 switch s
341 case -1 340 case -1
342 p_in = 2; 341 p_in = 2;
344 p_in = 3; 343 p_in = 3;
345 end 344 end
346 345
347 a = obj.gamma -1; 346 a = obj.gamma -1;
348 L = @(~,u,~)a*[0 -1/2*u 1]; 347 L = @(~,u,~)a*[0 -1/2*u 1];
349 g = @(t)[p_data(t)]; 348
350 349 closure = boundary_condition_L(obj, boundary, L, p_in);
351 350 end
352 closure = boundary_condition_L(obj, boundary, L, g, p_in); 351
353 352 % Return closure(q,[v; rho])
354 end 353 function closure = boundary_condition_inflow_rho(obj, boundary)
355
356 function closure = boundary_condition_inflow_rho(obj, boundary, rho_data, v_data)
357 [~,~,s] = obj.get_boundary_ops(boundary); 354 [~,~,s] = obj.get_boundary_ops(boundary);
358 355
359 switch s 356 switch s
360 case -1 357 case -1
361 p_in = [1 2]; 358 p_in = [1 2];
366 a = obj.gamma - 1; 363 a = obj.gamma - 1;
367 L = @(rho,~,~)[ 364 L = @(rho,~,~)[
368 0 1/rho 0; 365 0 1/rho 0;
369 1 0 0; 366 1 0 0;
370 ]; 367 ];
371 g = @(t)[ 368
372 v_data(t); 369 closure = boundary_condition_L(obj, boundary, L, p_in);
373 rho_data(t); 370 end
374 ]; 371
375 372 % Return closure(q,rho)
376 closure = boundary_condition_L(obj, boundary, L, g, p_in); 373 function closure = boundary_condition_outflow_rho(obj, boundary)
377 end
378
379 function closure = boundary_condition_outflow_rho(obj, boundary, rho_data)
380 [~,~,s] = obj.get_boundary_ops(boundary); 374 [~,~,s] = obj.get_boundary_ops(boundary);
381 375
382 switch s 376 switch s
383 case -1 377 case -1
384 p_in = 2; 378 p_in = 2;
385 case 1 379 case 1
386 p_in = 3; 380 p_in = 3;
387 end 381 end
388 382
389 L = @(~,~,~)[1 0 0]; 383 L = @(~,~,~)[1 0 0];
390 g = @(t)[rho_data(t)]; 384
391 385 closure = boundary_condition_L(obj, boundary, L, p_in);
392
393 closure = boundary_condition_L(obj, boundary, L, g, p_in);
394
395 end 386 end
396 387
397 % Set wall boundary condition v = 0. 388 % Set wall boundary condition v = 0.
398 function closure = boundary_condition_wall(obj,boundary) 389 function closure = boundary_condition_wall(obj,boundary)
399 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 390 [e_s,e_S,s] = obj.get_boundary_ops(boundary);