Mercurial > repos > public > sbplib
annotate +scheme/Euler1d.m @ 54:2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Wed, 11 Nov 2015 19:21:38 -0800 |
parents | 75ebf5d3cfe5 |
children | 9a647dcccbdd |
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 | 2 properties |
3 m % Number of points in each direction, possibly a vector | |
4 N % Number of points total | |
5 h % Grid spacing | |
6 u % Grid values | |
7 x % Values of x and y for each | |
8 order % Order accuracy for the approximation | |
9 | |
10 D % non-stabalized scheme operator | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
11 Fx |
0 | 12 M % Derivative norm |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
13 gamma |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
14 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
15 T |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
16 p |
50
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
17 c |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
18 |
0 | 19 |
20 H % Discrete norm | |
21 Hi | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
22 e_l, e_r, e_L, e_R; |
0 | 23 |
24 end | |
25 | |
26 methods | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
27 function obj = Euler1d(m,xlim,order,gama,opsGen,do_upwind) |
0 | 28 default_arg('opsGen',@sbp.Ordinary); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
29 default_arg('gama', 1.4); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
30 default_arg('do_upwind', false); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
31 gamma = gama; |
0 | 32 |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
33 [x, h] = util.get_grid(xlim{:},m); |
0 | 34 |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
35 if do_upwind |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
36 ops = sbp.Upwind(m,h,order); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
37 Dp = ops.derivatives.Dp; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
38 Dm = ops.derivatives.Dm; |
0 | 39 |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
40 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
|
41 Ddisp = (Dp - Dm)/2; |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
42 else |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
43 ops = opsGen(m,h,order); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
44 D1 = sparse(ops.derivatives.D1); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
45 end |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
46 |
0 | 47 H = sparse(ops.norms.H); |
48 Hi = sparse(ops.norms.HI); | |
49 e_l = sparse(ops.boundary.e_1); | |
50 e_r = sparse(ops.boundary.e_m); | |
51 | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
52 I_x = speye(m); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
53 I_3 = speye(3); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
54 |
0 | 55 D1 = kr(D1, I_3); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
56 if do_upwind |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
57 Ddisp = kr(Ddisp,I_3); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
58 end |
0 | 59 |
60 % Norms | |
61 obj.H = kr(H,I_3); | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
62 obj.Hi = kr(Hi,I_3); |
0 | 63 |
64 % Boundary operators | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
65 obj.e_l = e_l; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
66 obj.e_r = e_r; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
67 obj.e_L = kr(e_l,I_3); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
68 obj.e_R = kr(e_r,I_3); |
0 | 69 |
70 obj.m = m; | |
71 obj.h = h; | |
72 obj.order = order; | |
73 | |
74 % Man har Q_t+F_x=0 i 1D Euler, där | |
75 % q=[rho, rho*u, e]^T | |
76 % F=[rho*u, rho*u^2+p, (e+p)*u] ^T | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
77 % p=(gamma-1)*(e-rho*u^2/2); |
0 | 78 |
79 %Solving on form q_t + F_x = 0 | |
50
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
80 function o = F(Q) |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
81 % Flux: f = [q2; q2.^2/q1 + p(q); (q3+p(q))*q2/q1]; |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
82 o = [Q(2,:); Q(2,:).^2/Q(1,:) + p(Q); (Q(3,:)+p(Q)).*Q(2,:)./Q(1,:)]; |
0 | 83 end |
84 | |
85 % Equation of state | |
50
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
86 function o = p(Q) |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
87 % Pressure p = (gamma-1)*(q3-q2.^2/q1/2) |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
88 o = (gamma-1)*(Q(3,:)-Q(2,:).^2/Q(1,:)/2); |
0 | 89 end |
90 | |
50
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
91 function o = c(Q) |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
92 % Speed of light c = sqrt(obj.gamma*p/rho); |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
93 o = sqrt(gamma*p(Q)./Q(1,:)); |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
94 end |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
95 |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
96 function o = T(q) |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
97 % T is the transformation such that A = T*Lambda*inv(T) |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
98 % where Lambda=diag(u, u+c, u-c) |
0 | 99 rho = q(1); |
100 u = q(2)/q(1); | |
101 e = q(3); | |
102 | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
103 c = sqrt(gamma*p(q)/rho); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
104 |
0 | 105 sqrt2gamm = sqrt(2*(gamma-1)); |
106 | |
107 o = [ | |
108 sqrt2gamm*rho , rho , rho ; | |
109 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ; | |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
110 sqrt2gamm*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c , e+(gamma-1)*(e-rho*u^2/2)-rho*u*c ; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
111 ]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
112 % Devide columns by stuff to get rid of extra comp? |
0 | 113 end |
114 | |
115 function o = Fx(q) | |
50
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
116 Q = reshape(q,3,m); |
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
117 o = reshape(F(Q),3*m,1); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
118 o = D1*o; |
0 | 119 end |
120 | |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
121 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
|
122 Q = reshape(q,3,m); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
123 f = reshape(F(Q),3*m,1); |
0 | 124 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
125 c = c(Q); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
126 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
|
127 % 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
|
128 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
129 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
|
130 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
|
131 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
|
132 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
133 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
|
134 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
135 o = -D1*f + Ddisp*lamb_q; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
136 end |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
137 |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
138 if do_upwind |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
139 obj.D = @Fx_disp; |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
140 else |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
141 obj.D = @(q)-Fx(q); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
142 end |
0 | 143 |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
144 obj.Fx = @Fx; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
145 obj.T = @T; |
0 | 146 obj.u = x; |
147 obj.x = kr(x,ones(3,1)); | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
148 obj.p = @p; |
50
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
149 obj.c = @c; |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
150 obj.gamma = gamma; |
0 | 151 end |
152 | |
153 | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
154 % Enforces the boundary conditions |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
155 % w+ = R*w- + g(t) |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
156 function closure = boundary_condition(obj,boundary, type, varargin) |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
157 [e_s,e_S,s] = obj.get_boundary_ops(boundary); |
0 | 158 |
159 % Boundary condition on form | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
160 % w_in = R*w_out + g, where g is data |
0 | 161 |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
162 switch type |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
163 case 'wall' |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
164 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
|
165 case 'inflow' |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
166 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
|
167 case 'outflow' |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
168 closure = obj.boundary_condition_outflow(boundary,varargin{:}); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
169 otherwise |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
170 error('Unsupported bc type: %s', type); |
0 | 171 end |
172 | |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
173 end |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
174 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
175 function closure = boundary_condition_inflow(obj, boundary, p_data, v_data) |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
176 [e_s,e_S,s] = obj.get_boundary_ops(boundary); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
177 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
178 switch s |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
179 case -1 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
180 p_in = [1 2]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
181 p_ot = [3]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
182 case 1 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
183 p_in = [1 3]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
184 p_ot = [2]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
185 otherwise |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
186 error(); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
187 end |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
188 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
189 p = [p_in, p_ot]; % Permutation to sort |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
190 pt(p) = 1:length(p); % Inverse permutation |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
191 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
192 gamma = obj.gamma; |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
193 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
194 function o = closure_fun(q,t) |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
195 q_s = e_S'*q; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
196 rho = q_s(1); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
197 u = q_s(2)/rho; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
198 e = q_s(3); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
199 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
200 c = obj.c(q_s); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
201 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
202 T = obj.T(q_s); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
203 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
204 % L = [ |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
205 % 0 1 0; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
206 % 0 -1/2*u 1; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
207 % ]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
208 LTpi = [ |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
209 2^(1/2)/(2*rho*u*(gamma - 1)^(1/2)), -(2^(1/2)*(u - c*s))/(u*(2*e*gamma - rho*(gamma*u^2 - c*s*u + 2*c*s))*(gamma - 1)^(1/2)); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
210 0, 2/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2)); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
211 ]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
212 R = [ |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
213 (2^(1/2)*c*s*(rho*(gamma - 1)*u^2 + 2*rho*u - 2*e*gamma))/(u*(gamma - 1)^(1/2)*(- gamma*rho*u^2 + c*rho*s*u + 2*e*gamma - 2*c*rho*s)); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
214 (2*c*rho*s*(u - 2))/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2)) - 1; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
215 ]; |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
216 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
217 tau1 = -2*[ |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
218 u, 0; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
219 0, c + abs(u); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
220 ]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
221 tau2 = [0 0]; % Penalty only on ingoing char. |
0 | 222 |
223 | |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
224 tauHat = [tau1; tau2]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
225 tau = -s*e_S*sparse(T*tauHat(pt,:)); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
226 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
227 w_s = inv(T)*q_s; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
228 % w_s = T\q_s; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
229 % w_s = Tinv * q_s; % Med analytisk matris |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
230 w_in = w_s(p_in); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
231 w_ot = w_s(p_ot); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
232 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
233 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
234 o = obj.Hi * tau * (w_in - R*w_ot - LTpi*[v_data(t); p_data(t)]); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
235 end |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
236 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
237 closure = @closure_fun; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
238 end |
0 | 239 |
240 | |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
241 function closure = boundary_condition_outflow(obj, boundary, p_data) |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
242 [e_s,e_S,s] = obj.get_boundary_ops(boundary); |
0 | 243 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
244 switch s |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
245 case -1 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
246 p_in = 2; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
247 p_ot = [1 3]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
248 case 1 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
249 p_in = 3; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
250 p_ot = [1 2]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
251 otherwise |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
252 error(); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
253 end |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
254 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
255 p = [p_in, p_ot]; % Permutation to sort |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
256 pt(p) = 1:length(p); % Inverse permutation |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
257 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
258 gamma = obj.gamma; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
259 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
260 function o = closure_fun(q,t) |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
261 q_s = e_S'*q; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
262 rho = q_s(1); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
263 u = q_s(2)/rho; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
264 e = q_s(3); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
265 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
266 c = obj.c(q_s); |
0 | 267 |
268 | |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
269 T = obj.T(q_s); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
270 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
271 L = (gamma -1)*[0 -1/2*u 1]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
272 LTpi = 1/(L*T(:,p_in)); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
273 R = -[0 (gamma*(e - rho*u^2/2) - rho*u*c)/(gamma*(e - rho*u^2/2) + rho*u*c)]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
274 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
275 tau1 = [-2*(c - abs(u))]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
276 tau2 = [0; 0]; % Penalty only on ingoing char. |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
277 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
278 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
279 tauHat = [tau1; tau2]; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
280 tau = -s*e_S*sparse(T*tauHat(pt)); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
281 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
282 w_s = inv(T)*q_s; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
283 % w_s = T\q_s; |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
284 % w_s = Tinv * q_s; % Med analytisk matris |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
285 w_in = w_s(p_in); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
286 w_ot = w_s(p_ot); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
287 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
288 o = obj.Hi * tau * (w_in - R*w_ot - LTpi*p_data(t)); |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
289 end |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
290 |
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
291 closure = @closure_fun; |
0 | 292 |
293 end | |
294 | |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
295 % Set wall boundary condition v = 0. |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
296 function closure = boundary_condition_wall(obj,boundary) |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
297 [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
|
298 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
299 % v = 0 corresponds to |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
300 % L = [0 1 0]; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
301 % g = 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
302 % |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
303 % Tp = |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
304 % R = -(u-c)/(u+c) |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
305 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
306 % tau = alpha * (u+c) |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
307 % (alpha+1)(u+c) + 1/4* alpha^2|u-c| <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
308 % 4*(alpha+1)(u+c) + alpha^2|u-c| <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
309 % 4 * (alpha+1)(u+c) + alpha^2|u| + alpha^2*c <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
310 % |u|*(sgn(u)*4 + sgn(u)*4*alpha + alpha^2) + c*(4 + 4alpha + alpha^2) <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
311 % |u|*(alpha^2 + 4*sgn(u)*alpha + 4*sgn(u)) + c*(alpha+2)^2 <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
312 % |u|*[(alpha + 2*sgn(u))^2 - 4*(sgn(u)-1)] + c*(alpha+2)^2 <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
313 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
314 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
315 % om vi låtsas att u = 0: |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
316 % (alpha+1)c + 1/4 * alpha^2*c <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
317 % alpha^2 + 4*alpha +4 <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
318 % (alpha + 2)^2 <= 0 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
319 % alpha = -2 gives tau = -2*c; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
320 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
321 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
322 % 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
|
323 % vi låta den vara fri? |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
324 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
325 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
326 switch s |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
327 case -1 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
328 p_in = 2; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
329 p_zero = 1; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
330 p_ot = 3; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
331 case 1 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
332 p_in = 3; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
333 p_zero = 1; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
334 p_ot = 2; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
335 otherwise |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
336 error(); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
337 end |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
338 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
339 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
|
340 pt(p) = 1:length(p); % Inverse permutation |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
341 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
342 function o = closure_fun(q) |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
343 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
344 q_s = e_S'*q; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
345 rho = q_s(1); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
346 u = q_s(2)/rho; |
50
75ebf5d3cfe5
Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents:
49
diff
changeset
|
347 c = obj.c(q_s); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
348 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
349 T = obj.T(q_s); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
350 R = -(u-c)/(u+c); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
351 % l = [u, u+c, u-c]; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
352 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
353 % p_in = find(s*l <= 0); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
354 % p_ot = find(s*l > 0); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
355 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
356 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
357 tau1 = -2*c; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
358 tau2 = [0; 0]; % Penalty only on ingoing char. |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
359 |
54
2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents:
50
diff
changeset
|
360 % 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
|
361 % Lambda_ot = diag(abs(l(p_ot))); |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
362 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
363 tauHat = [tau1; tau2]; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
364 tau = -s*e_S*sparse(T*tauHat(pt)); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
365 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
366 w_s = inv(T)*q_s; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
367 % w_s = T\q_s; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
368 % w_s = Tinv * q_s; % Med analytisk matris |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
369 w_in = w_s(p_in); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
370 w_ot = w_s(p_ot); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
371 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
372 o = obj.Hi * tau * (w_in - R*w_ot); |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
373 end |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
374 |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
375 closure = @closure_fun; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
376 end |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
377 |
0 | 378 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
|
379 error('NOT DONE') |
0 | 380 % u denotes the solution in the own domain |
381 % v denotes the solution in the neighbour domain | |
382 [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); | |
383 [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
384 | |
385 tuning = 2; | |
386 | |
387 alpha_u = obj.alpha; | |
388 alpha_v = neighbour_scheme.alpha; | |
389 | |
390 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning; | |
391 % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning; | |
392 tau4 = s_u*alpha_u/2; | |
393 | |
394 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning; | |
395 sig3 = -s_u*alpha_u/2; | |
396 | |
397 phi2 = s_u*1/2; | |
398 | |
399 psi1 = -s_u*1/2; | |
400 | |
401 tau = tau1*e_u + tau4*d3_u; | |
402 sig = sig2*d1_u + sig3*d2_u ; | |
403 phi = phi2*d1_u ; | |
404 psi = psi1*e_u ; | |
405 | |
406 closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); | |
407 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); | |
408 end | |
409 | |
410 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
411 % 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
|
412 function [e,E,s] = get_boundary_ops(obj,boundary) |
0 | 413 switch boundary |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
414 case 'l' |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
415 e = obj.e_l; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
416 E = obj.e_L; |
0 | 417 s = -1; |
49
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
418 case 'r' |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
419 e = obj.e_r; |
8f0c2dc747dd
Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
420 E = obj.e_R; |
0 | 421 s = 1; |
422 otherwise | |
423 error('No such boundary: boundary = %s',boundary); | |
424 end | |
425 end | |
426 | |
427 function N = size(obj) | |
428 N = prod(obj.m); | |
429 end | |
430 | |
431 end | |
432 end |