Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 1062:e512714fb890 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/getBoundaryOp
author Martin Almquist <malmquist@stanford.edu>
date Mon, 14 Jan 2019 18:14:44 -0800
parents 2b1b944deae1
children 0c504a21432d
comparison
equal deleted inserted replaced
988:a72038b1f709 1062:e512714fb890
199 end 199 end
200 200
201 % Enforces the boundary conditions 201 % Enforces the boundary conditions
202 % w+ = R*w- + g(t) 202 % w+ = R*w- + g(t)
203 function closure = boundary_condition(obj,boundary, type, varargin) 203 function closure = boundary_condition(obj,boundary, type, varargin)
204 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 204 [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
205 s = obj.getBoundarySign(boundary);
205 206
206 % Boundary condition on form 207 % Boundary condition on form
207 % w_in = R*w_out + g, where g is data 208 % w_in = R*w_out + g, where g is data
208 209
209 switch type 210 switch type
230 % L = L(rho, u, e) 231 % L = L(rho, u, e)
231 % p_in are the indecies of the ingoing characteristics. 232 % p_in are the indecies of the ingoing characteristics.
232 % 233 %
233 % Returns closure(q,g) 234 % Returns closure(q,g)
234 function closure = boundary_condition_L(obj, boundary, L_fun, p_in) 235 function closure = boundary_condition_L(obj, boundary, L_fun, p_in)
235 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 236 [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
237 s = obj.getBoundarySign(boundary);
236 238
237 p_ot = 1:3; 239 p_ot = 1:3;
238 p_ot(p_in) = []; 240 p_ot(p_in) = [];
239 241
240 p = [p_in, p_ot]; % Permutation to sort 242 p = [p_in, p_ot]; % Permutation to sort
271 closure = @closure_fun; 273 closure = @closure_fun;
272 end 274 end
273 275
274 % Return closure(q,g) 276 % Return closure(q,g)
275 function closure = boundary_condition_char(obj,boundary) 277 function closure = boundary_condition_char(obj,boundary)
276 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 278 [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
279 s = obj.getBoundarySign(boundary);
277 280
278 function o = closure_fun(q, w_data) 281 function o = closure_fun(q, w_data)
279 q_s = e_S'*q; 282 q_s = e_S'*q;
280 rho = q_s(1); 283 rho = q_s(1);
281 u = q_s(2)/rho; 284 u = q_s(2)/rho;
312 end 315 end
313 316
314 317
315 % Return closure(q,[v; p]) 318 % Return closure(q,[v; p])
316 function closure = boundary_condition_inflow(obj, boundary) 319 function closure = boundary_condition_inflow(obj, boundary)
317 [~,~,s] = obj.get_boundary_ops(boundary); 320 s = obj.getBoundarySign(boundary);
318 321
319 switch s 322 switch s
320 case -1 323 case -1
321 p_in = [1 2]; 324 p_in = [1 2];
322 case 1 325 case 1
333 closure = @(q,p,v) closure_raw(q,[v; p]); 336 closure = @(q,p,v) closure_raw(q,[v; p]);
334 end 337 end
335 338
336 % Return closure(q, p) 339 % Return closure(q, p)
337 function closure = boundary_condition_outflow(obj, boundary) 340 function closure = boundary_condition_outflow(obj, boundary)
338 [~,~,s] = obj.get_boundary_ops(boundary); 341 s = obj.getBoundarySign(boundary);
339 342
340 switch s 343 switch s
341 case -1 344 case -1
342 p_in = 2; 345 p_in = 2;
343 case 1 346 case 1
350 closure = boundary_condition_L(obj, boundary, L, p_in); 353 closure = boundary_condition_L(obj, boundary, L, p_in);
351 end 354 end
352 355
353 % Return closure(q,[v; rho]) 356 % Return closure(q,[v; rho])
354 function closure = boundary_condition_inflow_rho(obj, boundary) 357 function closure = boundary_condition_inflow_rho(obj, boundary)
355 [~,~,s] = obj.get_boundary_ops(boundary); 358 s = obj.getBoundarySign(boundary);
356 359
357 switch s 360 switch s
358 case -1 361 case -1
359 p_in = [1 2]; 362 p_in = [1 2];
360 case 1 363 case 1
370 closure = boundary_condition_L(obj, boundary, L, p_in); 373 closure = boundary_condition_L(obj, boundary, L, p_in);
371 end 374 end
372 375
373 % Return closure(q,rho) 376 % Return closure(q,rho)
374 function closure = boundary_condition_outflow_rho(obj, boundary) 377 function closure = boundary_condition_outflow_rho(obj, boundary)
375 [~,~,s] = obj.get_boundary_ops(boundary); 378 s = obj.getBoundarySign(boundary);
376 379
377 switch s 380 switch s
378 case -1 381 case -1
379 p_in = 2; 382 p_in = 2;
380 case 1 383 case 1
386 closure = boundary_condition_L(obj, boundary, L, p_in); 389 closure = boundary_condition_L(obj, boundary, L, p_in);
387 end 390 end
388 391
389 % Set wall boundary condition v = 0. 392 % Set wall boundary condition v = 0.
390 function closure = boundary_condition_wall(obj,boundary) 393 function closure = boundary_condition_wall(obj,boundary)
391 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 394 [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary);
395 s = obj.getBoundarySign(boundary);
392 396
393 % Vill vi sätta penalty på karateristikan som är nära noll också eller vill 397 % Vill vi sätta penalty på karateristikan som är nära noll också eller vill
394 % vi låta den vara fri? 398 % vi låta den vara fri?
395 399
396 400
476 480
477 closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); 481 closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u');
478 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); 482 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
479 end 483 end
480 484
481 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 485 % Returns the boundary operator op for the boundary specified by the string boundary.
482 % The right boundary is considered the positive boundary 486 % op -- string or a cell array of strings
483 function [e,E,s] = get_boundary_ops(obj,boundary) 487 % boundary -- string
488 function varargout = getBoundaryOperator(obj, op, boundary)
489
490 if ~iscell(op)
491 op = {op};
492 end
493
494 for i = 1:numel(op)
495 switch op{i}
496 case 'e'
497 switch boundary
498 case 'l'
499 e = obj.e_l;
500 case 'r'
501 e = obj.e_r;
502 otherwise
503 error('No such boundary: boundary = %s',boundary);
504 end
505 varargout{i} = e;
506
507 case 'E'
508 switch boundary
509 case 'l'
510 E = obj.e_L;
511 case 'r'
512 E = obj.e_R;
513 otherwise
514 error('No such boundary: boundary = %s',boundary);
515 end
516 varargout{i} = E;
517 end
518 end
519 end
520
521 % Returns the boundary sign. The right boundary is considered the positive boundary
522 % boundary -- string
523 function s = getBoundarySign(obj, boundary)
484 switch boundary 524 switch boundary
485 case 'l' 525 case {'r'}
486 e = obj.e_l; 526 s = 1;
487 E = obj.e_L; 527 case {'l'}
488 s = -1; 528 s = -1;
489 case 'r'
490 e = obj.e_r;
491 E = obj.e_R;
492 s = 1;
493 otherwise 529 otherwise
494 error('No such boundary: boundary = %s',boundary); 530 error('No such boundary: boundary = %s',boundary);
495 end 531 end
496 end 532 end
497 533