Mercurial > repos > public > sbplib

annotate +scheme/Euler1d.m @ 1023:defc9d0cc1f2 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove incorrect assertion of the number of BC:s
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 07 Jan 2019 12:06:49 +0100
parents a35ed1d124d3
children 2b1b944deae1
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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=[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
80 % p=(gamma-1)*(e-rho*u^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 = -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
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/2*Q(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 = T*Lambda*inv(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 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ;
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
162 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 ;
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_S*sparse(T*tauHat(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/2*obj.Hi * tau * inv(L*Tin)*(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 = -s*e_S*sparse(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/2*obj.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/2*u*a 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/2*u 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 = -s*e_S*sparse(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/2*obj.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
944
a35ed1d124d3 Change from opts to type in a few schemes
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
449 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
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)/(2*delt_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 = tau1*e_u + tau4*d3_u;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
473 sig = sig2*d1_u + sig3*d2_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*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u');
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
478 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_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