sbplib: +scheme/Euler1d.m annotate

annotate +scheme/Euler1d.m @ 577:e45c9b56d50d feature/grids

Add an Empty grid class The need turned up for the flexural code when we may or may not have a grid for the open water and want to plot that solution. In case there is no open water we need an empty grid to plot the empty gridfunction against to avoid errors.

author	Jonatan Werpers <jonatan@werpers.com>
date	Thu, 07 Sep 2017 09:16:12 +0200
parents	f39f98b59f61
children	459eeb99130f

rev	line source
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	1 classdef Euler1d < scheme.Scheme
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	2 properties
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	3 m % Number of points in each direction, possibly a vector
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	4 N % Number of points total
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	5 h % Grid spacing
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	6 u % Grid values
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	7 x % Values of x and y for each
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	8 order % Order accuracy for the approximation
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	9
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	10 D % non-stabalized scheme operator
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	11 M % Derivative norm
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	12 gamma
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	13
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	14 H % Discrete norm
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	15 Hi
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	16 e_l, e_r, e_L, e_R;
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	17
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	18 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	19
93 5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	20 properties (Constant)
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	21 SUBSONIC_INFLOW = 1;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	22 SUBSONIC_OUTFLOW = -1;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	23 NO_FLOW = 0;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	24 SUPERSONIC_INFLOW = 2;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	25 SUPERSONIC_OUTFLOW = -2;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	26 end
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	27
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	28 methods
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	29 function obj = Euler1d(m,xlim,order,gama,opsGen,do_upwind)
302 f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	30 default_arg('opsGen',@sbp.D2Standard);
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	31 default_arg('gama', 1.4);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	32 default_arg('do_upwind', false);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	33 gamma = gama;
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	34
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	35 [x, h] = util.get_grid(xlim{:},m);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	36
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	37 if do_upwind
302 f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	38 ops = sbp.D1Upwind(m,xlim,order);
f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	39 Dp = ops.Dp;
f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	40 Dm = ops.Dm;
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	41
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	42 D1 = (Dp + Dm)/2;
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	43 Ddisp = (Dp - Dm)/2;
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	44 else
302 f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	45 ops = opsGen(m,xlim,order);
f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	46 printExpr('issparse(ops.D1)');
f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	47 D1 = ops.D1;
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	48 end
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	49
302 f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	50 H = sparse(ops.H);
f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	51 Hi = sparse(ops.HI);
f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	52 e_l = sparse(ops.e_l);
f39f98b59f61 Fixes in Upwind operators and Euler1D scheme. Jonatan Werpers <jonatan@werpers.com> parents: 111 diff changeset	53 e_r = sparse(ops.e_r);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	54
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	55 I_x = speye(m);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	56 I_3 = speye(3);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	57
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	58 D1 = kr(D1, I_3);
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	59 if do_upwind
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	60 Ddisp = kr(Ddisp,I_3);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	61 end
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	62
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	63 % Norms
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	64 obj.H = kr(H,I_3);
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	65 obj.Hi = kr(Hi,I_3);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	66
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	67 % Boundary operators
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	68 obj.e_l = e_l;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	69 obj.e_r = e_r;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	70 obj.e_L = kr(e_l,I_3);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	71 obj.e_R = kr(e_r,I_3);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	72
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	73 obj.m = m;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	74 obj.h = h;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	75 obj.order = order;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	76
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	77 % Man har Q_t+F_x=0 i 1D Euler, där
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	78 % q=[rho, rho*u, e]^T
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	79 % F=[rhou, rhou^2+p, (e+p)*u] ^T
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	80 % p=(gamma-1)(e-rhou^2/2);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	81
61 183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	82
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	83 %Solving on form q_t + F_x = 0
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	84
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	85 function o = Fx(q)
50 75ebf5d3cfe5 Performance improvement for Euler1d Jonatan Werpers <jonatan@werpers.com> parents: 49 diff changeset	86 Q = reshape(q,3,m);
61 183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	87 o = reshape(obj.F(Q),3*m,1);
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	88 o = D1*o;
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	89 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	90
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	91 function o = Fx_disp(q)
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	92 Q = reshape(q,3,m);
61 183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	93 f = reshape(obj.F(Q),3*m,1);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	94
61 183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	95 c = obj.c(Q);
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	96 lambda_max = c+abs(Q(2,:)./Q(1,:));
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	97 % lambda_max = max(lambda_max);
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	98
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	99 lamb_Q(1,:) = lambda_max.*Q(1,:);
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	100 lamb_Q(2,:) = lambda_max.*Q(2,:);
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	101 lamb_Q(3,:) = lambda_max.*Q(3,:);
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	102
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	103 lamb_q = reshape(lamb_Q,3*m, 1);
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	104
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	105 o = -D1f + Ddisplamb_q;
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	106 end
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	107
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	108 if do_upwind
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	109 obj.D = @Fx_disp;
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	110 else
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	111 obj.D = @(q)-Fx(q);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	112 end
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	113
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	114 obj.u = x;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	115 obj.x = kr(x,ones(3,1));
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	116 obj.gamma = gamma;
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	117 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	118
61 183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	119 % Flux function
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	120 function o = F(obj, Q)
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	121 % Flux: f = [q2; q2.^2/q1 + p(q); (q3+p(q))*q2/q1];
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	122 p = obj.p(Q);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	123 o = [Q(2,:); Q(2,:).^2./Q(1,:) + p; (Q(3,:)+p).*Q(2,:)./Q(1,:)];
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	124 end
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	125
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	126 % Equation of state
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	127 function o = p(obj, Q)
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	128 % Pressure p = (gamma-1)*(q3-q2.^2/q1/2)
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	129 o = (obj.gamma-1)( Q(3,:)-1/2Q(2,:).^2./Q(1,:) );
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	130 end
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	131
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	132 % Speed of sound
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	133 function o = c(obj, Q)
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	134 % Speed of light c = sqrt(obj.gamma*p/rho);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	135 o = sqrt(obj.gamma*obj.p(Q)./Q(1,:));
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	136 end
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	137
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	138 % Eigen value matrix
62 00344139cd7d Euler1d. Fixed mistake in last commit. Jonatan Werpers <jonatan@werpers.com> parents: 61 diff changeset	139 function o = Lambda(obj, q)
61 183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	140 u = q(2)/q(1);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	141 c = obj.c(q);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	142 L = [u, u+c, u-c];
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	143 o = diag(L);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	144 end
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	145
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	146 % Diagonalization transformation
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	147 function o = T(obj, q)
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	148 % T is the transformation such that A = TLambdainv(T)
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	149 % where Lambda=diag(u, u+c, u-c)
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	150 rho = q(1);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	151 u = q(2)/q(1);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	152 e = q(3);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	153 gamma = obj.gamma;
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	154
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	155 c = sqrt(gamma*obj.p(q)/rho);
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	156
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	157 sqrt2gamm = sqrt(2*(gamma-1));
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	158
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	159 o = [
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	160 sqrt2gamm*rho , rho , rho ;
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	161 sqrt2gammrhou , rho(u+c) , rho(u-c) ;
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	162 sqrt2gammrhou^2/2, e+(gamma-1)(e-rhou^2/2)+rhouc , e+(gamma-1)(e-rhou^2/2)-rhouc ;
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	163 ];
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	164 % Devide columns by stuff to get rid of extra comp?
183d9349b4c1 Refactored euler1d to have real methods. Jonatan Werpers <jonatan@werpers.com> parents: 60 diff changeset	165 end
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	166
93 5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	167 function fs = flowStateL(obj, q)
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	168 q_l = obj.e_L'*q;
111 0e66299592cc Fixed some bugs regarding flowState determination. Jonatan Werpers <jonatan@werpers.com> parents: 110 diff changeset	169 c = obj.c(q_l);
93 5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	170 v = q_l(2,:)/q_l(1,:);
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	171
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	172 if v > c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	173 fs = scheme.Euler1d.SUPERSONIC_INFLOW;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	174 elseif v > 0
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	175 fs = scheme.Euler1d.SUBSONIC_INFLOW;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	176 elseif v > -c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	177 fs = scheme.Euler1d.SUBSONIC_OUTFLOW;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	178 else
111 0e66299592cc Fixed some bugs regarding flowState determination. Jonatan Werpers <jonatan@werpers.com> parents: 110 diff changeset	179 fs = scheme.Euler1d.SUPERSONIC_OUTFLOW;
93 5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	180 end
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	181 end
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	182
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	183 % returns positiv values for inlfow, negative for outflow.
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	184 % +-1 for subsonic
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	185 function fs = flowStateR(obj, q)
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	186 q_r = obj.e_R'*q;
111 0e66299592cc Fixed some bugs regarding flowState determination. Jonatan Werpers <jonatan@werpers.com> parents: 110 diff changeset	187 c = obj.c(q_r);
93 5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	188 v = q_r(2,:)/q_r(1,:);
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	189
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	190 if v < -c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	191 fs = scheme.Euler1d.SUPERSONIC_INFLOW;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	192 elseif v < 0
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	193 fs = scheme.Euler1d.SUBSONIC_INFLOW;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	194 elseif v < c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	195 fs = scheme.Euler1d.SUBSONIC_OUTFLOW;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	196 else
111 0e66299592cc Fixed some bugs regarding flowState determination. Jonatan Werpers <jonatan@werpers.com> parents: 110 diff changeset	197 fs = scheme.Euler1d.SUPERSONIC_OUTFLOW;
93 5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	198 end
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	199 end
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary. Jonatan Werpers <jonatan@werpers.com> parents: 92 diff changeset	200
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	201 % Enforces the boundary conditions
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	202 % w+ = R*w- + g(t)
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	203 function closure = boundary_condition(obj,boundary, type, varargin)
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	204 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	205
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	206 % Boundary condition on form
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	207 % w_in = R*w_out + g, where g is data
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	208
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	209 switch type
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	210 case 'wall'
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	211 closure = obj.boundary_condition_wall(boundary);
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	212 case 'inflow'
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	213 closure = obj.boundary_condition_inflow(boundary,varargin{:});
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	214 case 'outflow'
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	215 closure = obj.boundary_condition_outflow(boundary,varargin{:});
91 2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	216 case 'inflow_rho'
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	217 closure = obj.boundary_condition_inflow_rho(boundary,varargin{:});
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	218 case 'outflow_rho'
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	219 closure = obj.boundary_condition_outflow_rho(boundary,varargin{:});
57 9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	220 case 'char'
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	221 closure = obj.boundary_condition_char(boundary,varargin{:});
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	222 otherwise
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	223 error('Unsupported bc type: %s', type);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	224 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	225
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	226 end
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	227
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	228
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	229 % Sets the boundary condition Lq = g, where
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	230 % L = L(rho, u, e)
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	231 % p_in are the indecies of the ingoing characteristics.
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	232 %
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	233 % Returns closure(q,g)
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	234 function closure = boundary_condition_L(obj, boundary, L_fun, p_in)
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	235 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	236
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	237 p_ot = 1:3;
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	238 p_ot(p_in) = [];
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	239
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	240 p = [p_in, p_ot]; % Permutation to sort
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	241 pt(p) = 1:length(p); % Inverse permutation
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	242
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	243 function o = closure_fun(q,g)
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	244 % Extract solution at the boundary
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	245 q_s = e_S'*q;
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	246 rho = q_s(1);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	247 u = q_s(2)/rho;
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	248 e = q_s(3);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	249
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	250 c = obj.c(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	251
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	252 % Calculate transformation matrix
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	253 T = obj.T(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	254 Tin = T(:,p_in);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	255 Tot = T(:,p_ot);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	256
92 ed7c7d651428 Euler1d: Fixed sign bug. Jonatan Werpers <jonatan@werpers.com> parents: 91 diff changeset	257 % Calculate eigen value matrix
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	258 Lambda = obj.Lambda(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	259
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	260 % Setup the penalty parameter
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	261 tau1 = -2*abs(Lambda(p_in,p_in));
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	262 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	263
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	264 tauHat = [tau1; tau2];
92 ed7c7d651428 Euler1d: Fixed sign bug. Jonatan Werpers <jonatan@werpers.com> parents: 91 diff changeset	265 tau = e_Ssparse(TtauHat(pt,:));
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	266
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	267 L = L_fun(rho,u,e);
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	268
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	269 o = 1/2obj.Hi tau * inv(LTin)(L*q_s - g);
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	270 end
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	271 closure = @closure_fun;
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	272 end
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	273
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	274 % Return closure(q,g)
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	275 function closure = boundary_condition_char(obj,boundary)
57 9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	276 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	277
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	278 function o = closure_fun(q, w_data)
57 9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	279 q_s = e_S'*q;
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	280 rho = q_s(1);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	281 u = q_s(2)/rho;
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	282 e = q_s(3);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	283
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	284 c = obj.c(q_s);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	285
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	286 Lambda = [u, u+c, u-c];
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	287
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	288 p_in = find(s*Lambda < 0);
100 ce4eecbcb915 Made char BC in euler more robust. Jonatan Werpers <jonatan@werpers.com> parents: 93 diff changeset	289 p_ot = 1:3;
ce4eecbcb915 Made char BC in euler more robust. Jonatan Werpers <jonatan@werpers.com> parents: 93 diff changeset	290 p_ot(p_in) = [];
57 9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	291 p = [p_in p_ot];
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	292 pt(p) = 1:length(p);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	293
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	294 T = obj.T(q_s);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	295
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	296 tau1 = -2*diag(abs(Lambda(p_in)));
57 9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	297 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	298
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	299 tauHat = [tau1; tau2];
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	300
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	301 tau = -se_Ssparse(T*tauHat(pt,:));
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	302
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	303 w_s = inv(T)*q_s;
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	304 w_in = w_s(p_in);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	305
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	306 w_in_data = w_data(p_in);
57 9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	307
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	308 o = 1/2obj.Hi tau * (w_in - w_in_data);
57 9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	309 end
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	310
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	311 closure = @closure_fun;
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	312 end
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 54 diff changeset	313
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	314
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	315 % Return closure(q,[v; p])
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	316 function closure = boundary_condition_inflow(obj, boundary)
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	317 [~,~,s] = obj.get_boundary_ops(boundary);
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	318
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	319 switch s
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	320 case -1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	321 p_in = [1 2];
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	322 case 1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	323 p_in = [1 3];
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	324 end
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	325
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	326 a = obj.gamma - 1;
63 21b18edb9667 Fixed error in bc. Jonatan Werpers <jonatan@werpers.com> parents: 62 diff changeset	327 L = @(rho,u,~)[
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	328 0 1/rho 0; %v
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	329 0 -1/2ua a; %p
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	330 ];
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	331
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	332 closure_raw = boundary_condition_L(obj, boundary, L, g, p_in);
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	333 closure = @(q,p,v) closure_raw(q,[v; p]);
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	334 end
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	335
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	336 % Return closure(q, p)
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	337 function closure = boundary_condition_outflow(obj, boundary)
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	338 [~,~,s] = obj.get_boundary_ops(boundary);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	339
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	340 switch s
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	341 case -1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	342 p_in = 2;
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	343 case 1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	344 p_in = 3;
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	345 end
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	346
60 e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	347 a = obj.gamma -1;
e707cd9e6d0a Euler1d: Created a general bc routine. Jonatan Werpers <jonatan@werpers.com> parents: 59 diff changeset	348 L = @(~,u,~)a[0 -1/2u 1];
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	349
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	350 closure = boundary_condition_L(obj, boundary, L, p_in);
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	351 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	352
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	353 % Return closure(q,[v; rho])
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	354 function closure = boundary_condition_inflow_rho(obj, boundary)
91 2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	355 [~,~,s] = obj.get_boundary_ops(boundary);
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	356
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	357 switch s
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	358 case -1
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	359 p_in = [1 2];
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	360 case 1
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	361 p_in = [1 3];
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	362 end
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	363
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	364 a = obj.gamma - 1;
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	365 L = @(rho,~,~)[
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	366 0 1/rho 0;
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	367 1 0 0;
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	368 ];
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	369
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	370 closure = boundary_condition_L(obj, boundary, L, p_in);
91 2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	371 end
2102af217134 Added inflow bc for rho Jonatan Werpers <jonatan@werpers.com> parents: 78 diff changeset	372
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	373 % Return closure(q,rho)
f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	374 function closure = boundary_condition_outflow_rho(obj, boundary)
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	375 [~,~,s] = obj.get_boundary_ops(boundary);
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	376
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	377 switch s
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	378 case -1
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	379 p_in = 2;
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	380 case 1
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	381 p_in = 3;
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	382 end
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	383
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	384 L = @(~,~,~)[1 0 0];
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	385
110 f5ed7ff58115 Changed so that closures accepts bc data instead of the closure cretaorChanged so that closures accepts bc data instead of the closure cretaors. Jonatan Werpers <jonatan@werpers.com> parents: 100 diff changeset	386 closure = boundary_condition_L(obj, boundary, L, p_in);
77 0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	387 end
0b07ff8a0a12 Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2. Jonatan Werpers <jonatan@werpers.com> parents: 63 diff changeset	388
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	389 % Set wall boundary condition v = 0.
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	390 function closure = boundary_condition_wall(obj,boundary)
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	391 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	392
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	393 % Vill vi sätta penalty på karateristikan som är nära noll också eller vill
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	394 % vi låta den vara fri?
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	395
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	396
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	397 switch s
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	398 case -1
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	399 p_in = 2;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	400 p_zero = 1;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	401 p_ot = 3;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	402 case 1
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	403 p_in = 3;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	404 p_zero = 1;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	405 p_ot = 2;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	406 otherwise
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	407 error();
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	408 end
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	409
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	410 p = [p_in, p_zero, p_ot]; % Permutation to sort
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	411 pt(p) = 1:length(p); % Inverse permutation
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	412
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	413 function o = closure_fun(q)
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	414
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	415 q_s = e_S'*q;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	416 rho = q_s(1);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	417 u = q_s(2)/rho;
50 75ebf5d3cfe5 Performance improvement for Euler1d Jonatan Werpers <jonatan@werpers.com> parents: 49 diff changeset	418 c = obj.c(q_s);
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	419
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	420 T = obj.T(q_s);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	421 R = -(u-c)/(u+c);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	422 % l = [u, u+c, u-c];
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	423
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	424 % p_in = find(s*l <= 0);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	425 % p_ot = find(s*l > 0);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	426
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	427
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	428 tau1 = -2*c;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	429 tau2 = [0; 0]; % Penalty only on ingoing char.
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	430
54 2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	431 % Lambda_in = diag(abs(l(p_in)));
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC. Jonatan Werpers <jonatan@werpers.com> parents: 50 diff changeset	432 % Lambda_ot = diag(abs(l(p_ot)));
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	433
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	434 tauHat = [tau1; tau2];
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	435 tau = -se_Ssparse(T*tauHat(pt));
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	436
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	437 w_s = inv(T)*q_s;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	438 % w_s = T\q_s;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	439 % w_s = Tinv * q_s; % Med analytisk matris
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	440 w_in = w_s(p_in);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	441 w_ot = w_s(p_ot);
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	442
78 80948a4084f3 Added the factor or 1/2 to the wall bc term as well. Jonatan Werpers <jonatan@werpers.com> parents: 77 diff changeset	443 o = 1/2obj.Hi tau * (w_in - R*w_ot);
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	444 end
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	445
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	446 closure = @closure_fun;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	447 end
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	448
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	449 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	450 error('NOT DONE')
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	451 % u denotes the solution in the own domain
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	452 % v denotes the solution in the neighbour domain
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	453 [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	454 [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	455
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	456 tuning = 2;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	457
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	458 alpha_u = obj.alpha;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	459 alpha_v = neighbour_scheme.alpha;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	460
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	461 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	462 % tau1 = (alpha_u/2 + alpha_v/2)/(2delt_u)tuning;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	463 tau4 = s_u*alpha_u/2;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	464
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	465 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	466 sig3 = -s_u*alpha_u/2;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	467
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	468 phi2 = s_u*1/2;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	469
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	470 psi1 = -s_u*1/2;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	471
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	472 tau = tau1e_u + tau4d3_u;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	473 sig = sig2d1_u + sig3d2_u ;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	474 phi = phi2*d1_u ;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	475 psi = psi1*e_u ;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	476
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	477 closure = halfnorm_inv(taue_u' + sigd1_u' + phialpha_ud2_u' + psialpha_u*d3_u');
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	478 penalty = -halfnorm_inv(taue_v' + sigd1_v' + phialpha_vd2_v' + psialpha_v*d3_v');
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	479 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	480
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	481 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	482 % The right boundary is considered the positive boundary
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	483 function [e,E,s] = get_boundary_ops(obj,boundary)
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	484 switch boundary
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	485 case 'l'
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	486 e = obj.e_l;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	487 E = obj.e_L;
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	488 s = -1;
49 8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	489 case 'r'
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	490 e = obj.e_r;
8f0c2dc747dd Made lots of updates to Euler1d. Jonatan Werpers <jonatan@werpers.com> parents: 0 diff changeset	491 E = obj.e_R;
0 48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	492 s = 1;
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	493 otherwise
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	494 error('No such boundary: boundary = %s',boundary);
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	495 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	496 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	497
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	498 function N = size(obj)
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	499 N = prod(obj.m);
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	500 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	501
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	502 end
48b6fb693025 Initial commit. Jonatan Werpers <jonatan@werpers.com> parents: diff changeset	503 end

Mercurial > repos > public > sbplib

annotate +scheme/Euler1d.m @ 577:e45c9b56d50d feature/grids