Mercurial > repos > public > sbplib

annotate +scheme/Euler1d.m @ 100:ce4eecbcb915

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Made char BC in euler more robust.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 04 Dec 2015 16:35:30 +0100
parents 5c41941b9e5e
children f5ed7ff58115
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)
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 default_arg('opsGen',@sbp.Ordinary);
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
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
38 ops = sbp.Upwind(m,h,order);
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
39 Dp = ops.derivatives.Dp;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
40 Dm = ops.derivatives.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
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
45 ops = opsGen(m,h,order);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
46 D1 = sparse(ops.derivatives.D1);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
47 end
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
48
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 H = sparse(ops.norms.H);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 Hi = sparse(ops.norms.HI);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 e_l = sparse(ops.boundary.e_1);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 e_r = sparse(ops.boundary.e_m);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
54 I_x = speye(m);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
55 I_3 = speye(3);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
56
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 D1 = kr(D1, I_3);
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
58 if do_upwind
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
59 Ddisp = kr(Ddisp,I_3);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
60 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 % Norms
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 obj.H = kr(H,I_3);
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
64 obj.Hi = kr(Hi,I_3);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 % Boundary operators
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
67 obj.e_l = e_l;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
68 obj.e_r = e_r;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
69 obj.e_L = kr(e_l,I_3);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
70 obj.e_R = kr(e_r,I_3);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 obj.m = m;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73 obj.h = h;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 obj.order = order;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 % Man har Q_t+F_x=0 i 1D Euler, där
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 % q=[rho, rho*u, e]^T
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 % 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
79 % p=(gamma-1)*(e-rho*u^2/2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80
61
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
81
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 %Solving on form q_t + F_x = 0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84 function o = Fx(q)
50
75ebf5d3cfe5 Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents: 49
diff changeset
85 Q = reshape(q,3,m);
61
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
86 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
87 o = D1*o;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
90 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
91 Q = reshape(q,3,m);
61
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
92 f = reshape(obj.F(Q),3*m,1);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93
61
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
94 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
95 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
96 % 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
97
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
98 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
99 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
100 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
101
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
102 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
103
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
104 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
105 end
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
106
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
107 if do_upwind
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
108 obj.D = @Fx_disp;
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
109 else
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
110 obj.D = @(q)-Fx(q);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
111 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
113 obj.u = x;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114 obj.x = kr(x,ones(3,1));
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
115 obj.gamma = gamma;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
116 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
117
61
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
118 % Flux function
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
119 function o = F(obj, Q)
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
120 % 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
121 p = obj.p(Q);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
122 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
123 end
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
124
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
125 % Equation of state
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
126 function o = p(obj, Q)
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
127 % 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
128 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
129 end
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
130
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
131 % Speed of sound
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
132 function o = c(obj, Q)
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
133 % 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
134 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
135 end
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
136
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
137 % Eigen value matrix
62
00344139cd7d Euler1d. Fixed mistake in last commit.
Jonatan Werpers <jonatan@werpers.com>
parents: 61
diff changeset
138 function o = Lambda(obj, q)
61
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
139 u = q(2)/q(1);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
140 c = obj.c(q);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
141 L = [u, u+c, u-c];
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
142 o = diag(L);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
143 end
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
144
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
145 % Diagonalization transformation
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
146 function o = T(obj, q)
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
147 % 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
148 % where Lambda=diag(u, u+c, u-c)
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
149 rho = q(1);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
150 u = q(2)/q(1);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
151 e = q(3);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
152 gamma = obj.gamma;
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
153
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
154 c = sqrt(gamma*obj.p(q)/rho);
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
155
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
156 sqrt2gamm = sqrt(2*(gamma-1));
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
157
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
158 o = [
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
159 sqrt2gamm*rho , rho , rho ;
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
160 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
161 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
162 ];
183d9349b4c1 Refactored euler1d to have real methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
163 % 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
164 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
165
93
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
166 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
167 q_l = obj.e_L'*q;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
168 c = obj.c(q);
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
169 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
170
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
171 if v > c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
172 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
173 elseif v > 0
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
174 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
175 elseif v > -c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
176 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
177 else
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
178 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
179 end
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
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
182 % 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
183 % +-1 for subsonic
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
184 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
185 q_r = obj.e_R'*q;
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
186 c = obj.c(q);
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
187 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
188
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
189 if v < -c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
190 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
191 elseif v < 0
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
192 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
193 elseif v < c
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
194 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
195 else
5c41941b9e5e Euler1D: added function to detect type of flow at a baoundary.
Jonatan Werpers <jonatan@werpers.com>
parents: 92
diff changeset
196 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
197 end
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
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
200 % Enforces the boundary conditions
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
201 % w+ = R*w- + g(t)
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
202 function closure = boundary_condition(obj,boundary, type, varargin)
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
203 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
204
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
205 % Boundary condition on form
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
206 % w_in = R*w_out + g, where g is data
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
207
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
208 switch type
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
209 case 'wall'
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
210 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
211 case 'inflow'
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
212 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
213 case 'outflow'
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
214 closure = obj.boundary_condition_outflow(boundary,varargin{:});
91
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
215 case 'inflow_rho'
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
216 closure = obj.boundary_condition_inflow_rho(boundary,varargin{:});
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
217 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
218 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
219 case 'char'
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
220 closure = obj.boundary_condition_char(boundary,varargin{:});
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
221 otherwise
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
222 error('Unsupported bc type: %s', type);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
223 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
224
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
225 end
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
226
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
227
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
228 % Sets the boundary condition Lq = g, where
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
229 % L = L(rho, u, e), g = g(t)
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
230 % p_in are the indecies of the ingoing characteristics.
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
231 function closure = boundary_condition_L(obj, boundary, L_fun, g_fun, p_in)
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
232 [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
233
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
234 p_ot = 1:3;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
235 p_ot(p_in) = [];
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 = [p_in, p_ot]; % Permutation to sort
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
238 pt(p) = 1:length(p); % Inverse permutation
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 function o = closure_fun(q,t)
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
241 % Extract solution at the boundary
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
242 q_s = e_S'*q;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
243 rho = q_s(1);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
244 u = q_s(2)/rho;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
245 e = q_s(3);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
246
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
247 c = obj.c(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
248
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
249 % Calculate transformation matrix
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
250 T = obj.T(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
251 Tin = T(:,p_in);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
252 Tot = T(:,p_ot);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
253
92
ed7c7d651428 Euler1d: Fixed sign bug.
Jonatan Werpers <jonatan@werpers.com>
parents: 91
diff changeset
254 % Calculate eigen value matrix
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
255 Lambda = obj.Lambda(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
256
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
257 % 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
258 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
259 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
260
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
261 tauHat = [tau1; tau2];
92
ed7c7d651428 Euler1d: Fixed sign bug.
Jonatan Werpers <jonatan@werpers.com>
parents: 91
diff changeset
262 tau = e_S*sparse(T*tauHat(pt,:));
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
263
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
264 L = L_fun(rho,u,e);
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
265 g = g_fun(t);
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 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
268 end
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
269 closure = @closure_fun;
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
57
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
272 function closure = boundary_condition_char(obj,boundary,w_data)
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
273 [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
274
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
275 function o = closure_fun(q,t)
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
276 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
277 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
278 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
279 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
280
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
281 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
282
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
283 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
284
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
285 p_in = find(s*Lambda < 0);
100
ce4eecbcb915 Made char BC in euler more robust.
Jonatan Werpers <jonatan@werpers.com>
parents: 93
diff changeset
286 p_ot = 1:3;
ce4eecbcb915 Made char BC in euler more robust.
Jonatan Werpers <jonatan@werpers.com>
parents: 93
diff changeset
287 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
288 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
289 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
290
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
291 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
292
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
293 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
294 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
295
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
296 tauHat = [tau1; tau2];
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
297
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
298 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
299
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
300 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
301 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
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_data = w_data(t);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
304 w_in_data = w_s_data(p_in);
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
305
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
306 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
307 end
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
308
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
309 closure = @closure_fun;
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 end
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
312
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
313 function closure = boundary_condition_inflow(obj, boundary, p_data, v_data)
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
314 [~,~,s] = obj.get_boundary_ops(boundary);
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
315
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
316 switch s
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
317 case -1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
318 p_in = [1 2];
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
319 case 1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
320 p_in = [1 3];
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
321 end
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
322
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
323 a = obj.gamma - 1;
63
21b18edb9667 Fixed error in bc.
Jonatan Werpers <jonatan@werpers.com>
parents: 62
diff changeset
324 L = @(rho,u,~)[
21b18edb9667 Fixed error in bc.
Jonatan Werpers <jonatan@werpers.com>
parents: 62
diff changeset
325 0 1/rho 0;
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
326 0 -1/2*u*a a;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
327 ];
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
328 g = @(t)[
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
329 v_data(t);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
330 p_data(t);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
331 ];
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
332
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
333 closure = boundary_condition_L(obj, boundary, L, g, p_in);
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
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
336
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
337 function closure = boundary_condition_outflow(obj, boundary, p_data)
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];
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
349 g = @(t)[p_data(t)];
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
350
59
e431c1260f52 Euler1d: Changed penalties to be calcualted numerically.
Jonatan Werpers <jonatan@werpers.com>
parents: 58
diff changeset
351
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
352 closure = boundary_condition_L(obj, boundary, L, g, p_in);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
353
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
354 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
355
91
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
356 function closure = boundary_condition_inflow_rho(obj, boundary, rho_data, v_data)
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
357 [~,~,s] = obj.get_boundary_ops(boundary);
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
358
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
359 switch s
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 2];
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
362 case 1
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
363 p_in = [1 3];
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
364 end
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
365
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
366 a = obj.gamma - 1;
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
367 L = @(rho,~,~)[
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
368 0 1/rho 0;
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
369 1 0 0;
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
370 ];
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
371 g = @(t)[
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
372 v_data(t);
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
373 rho_data(t);
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
374 ];
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
375
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
376 closure = boundary_condition_L(obj, boundary, L, g, p_in);
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
377 end
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
378
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
379 function closure = boundary_condition_outflow_rho(obj, boundary, rho_data)
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 [~,~,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
381
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 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
383 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
384 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
385 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
386 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
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
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
389 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
390 g = @(t)[rho_data(t)];
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
391
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
392
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
393 closure = boundary_condition_L(obj, boundary, L, g, p_in);
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
394
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
395 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
396
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
397 % Set wall boundary condition v = 0.
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
398 function closure = boundary_condition_wall(obj,boundary)
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
399 [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
400
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
401 % 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
402 % vi låta den vara fri?
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
403
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
404
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
405 switch s
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
406 case -1
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
407 p_in = 2;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
408 p_zero = 1;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
409 p_ot = 3;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
410 case 1
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
411 p_in = 3;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
412 p_zero = 1;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
413 p_ot = 2;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
414 otherwise
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
415 error();
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
416 end
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
417
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
418 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
419 pt(p) = 1:length(p); % Inverse permutation
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
420
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
421 function o = closure_fun(q)
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
422
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
423 q_s = e_S'*q;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
424 rho = q_s(1);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
425 u = q_s(2)/rho;
50
75ebf5d3cfe5 Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents: 49
diff changeset
426 c = obj.c(q_s);
49
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 T = obj.T(q_s);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
429 R = -(u-c)/(u+c);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
430 % l = [u, u+c, u-c];
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
431
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
432 % p_in = find(s*l <= 0);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
433 % p_ot = find(s*l > 0);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
434
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
435
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
436 tau1 = -2*c;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
437 tau2 = [0; 0]; % Penalty only on ingoing char.
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
438
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
439 % 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
440 % Lambda_ot = diag(abs(l(p_ot)));
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
441
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
442 tauHat = [tau1; tau2];
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
443 tau = -s*e_S*sparse(T*tauHat(pt));
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
444
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
445 w_s = inv(T)*q_s;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
446 % w_s = T\q_s;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
447 % w_s = Tinv * q_s; % Med analytisk matris
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
448 w_in = w_s(p_in);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
449 w_ot = w_s(p_ot);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
450
78
80948a4084f3 Added the factor or 1/2 to the wall bc term as well.
Jonatan Werpers <jonatan@werpers.com>
parents: 77
diff changeset
451 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
452 end
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
453
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
454 closure = @closure_fun;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
455 end
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
456
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
457 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
458 error('NOT DONE')
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
459 % u denotes the solution in the own domain
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
460 % v denotes the solution in the neighbour domain
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
461 [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
462 [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
463
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
464 tuning = 2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
465
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
466 alpha_u = obj.alpha;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
467 alpha_v = neighbour_scheme.alpha;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
468
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
469 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
470 % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
471 tau4 = s_u*alpha_u/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
472
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
473 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
474 sig3 = -s_u*alpha_u/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
475
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
476 phi2 = s_u*1/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
477
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
478 psi1 = -s_u*1/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
479
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
480 tau = tau1*e_u + tau4*d3_u;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
481 sig = sig2*d1_u + sig3*d2_u ;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
482 phi = phi2*d1_u ;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
483 psi = psi1*e_u ;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
484
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
485 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
486 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
487 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
488
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
489 % 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
490 % 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
491 function [e,E,s] = get_boundary_ops(obj,boundary)
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
492 switch boundary
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
493 case 'l'
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
494 e = obj.e_l;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
495 E = obj.e_L;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
496 s = -1;
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
497 case 'r'
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
498 e = obj.e_r;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
499 E = obj.e_R;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
500 s = 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
501 otherwise
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
502 error('No such boundary: boundary = %s',boundary);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
503 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
504 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
505
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
506 function N = size(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
507 N = prod(obj.m);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
508 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
509
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
510 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
511 end