Mercurial > repos > public > sbplib

annotate +scheme/Euler1d.m @ 311:713b125038a3 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed function names in renamed files.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 14:59:55 +0200
parents 0e66299592cc
children f39f98b59f61
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;
111
0e66299592cc Fixed some bugs regarding flowState determination.
Jonatan Werpers <jonatan@werpers.com>
parents: 110
diff changeset
168 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
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
111
0e66299592cc Fixed some bugs regarding flowState determination.
Jonatan Werpers <jonatan@werpers.com>
parents: 110
diff changeset
178 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
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;
111
0e66299592cc Fixed some bugs regarding flowState determination.
Jonatan Werpers <jonatan@werpers.com>
parents: 110
diff changeset
186 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
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
111
0e66299592cc Fixed some bugs regarding flowState determination.
Jonatan Werpers <jonatan@werpers.com>
parents: 110
diff changeset
196 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
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
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
229 % L = L(rho, u, e)
60
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.
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
231 %
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 % 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
233 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
234 [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
235
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
236 p_ot = 1:3;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
237 p_ot(p_in) = [];
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
238
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
239 p = [p_in, p_ot]; % Permutation to sort
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
240 pt(p) = 1:length(p); % Inverse permutation
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
241
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
242 function o = closure_fun(q,g)
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
243 % Extract solution at the boundary
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
244 q_s = e_S'*q;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
245 rho = q_s(1);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
246 u = q_s(2)/rho;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
247 e = q_s(3);
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 c = obj.c(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
250
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
251 % Calculate transformation matrix
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
252 T = obj.T(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
253 Tin = T(:,p_in);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
254 Tot = T(:,p_ot);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
255
92
ed7c7d651428 Euler1d: Fixed sign bug.
Jonatan Werpers <jonatan@werpers.com>
parents: 91
diff changeset
256 % Calculate eigen value matrix
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
257 Lambda = obj.Lambda(q_s);
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
258
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
259 % 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
260 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
261 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
262
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
263 tauHat = [tau1; tau2];
92
ed7c7d651428 Euler1d: Fixed sign bug.
Jonatan Werpers <jonatan@werpers.com>
parents: 91
diff changeset
264 tau = e_S*sparse(T*tauHat(pt,:));
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
265
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
266 L = L_fun(rho,u,e);
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
267
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
268 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
269 end
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
270 closure = @closure_fun;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
271 end
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
272
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
273 % 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
274 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
275 [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
276
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
277 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
278 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
279 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
280 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
281 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
282
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
283 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
284
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
285 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
286
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
287 p_in = find(s*Lambda < 0);
100
ce4eecbcb915 Made char BC in euler more robust.
Jonatan Werpers <jonatan@werpers.com>
parents: 93
diff changeset
288 p_ot = 1:3;
ce4eecbcb915 Made char BC in euler more robust.
Jonatan Werpers <jonatan@werpers.com>
parents: 93
diff changeset
289 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
290 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
291 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
292
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
293 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
294
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
295 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
296 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
297
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
298 tauHat = [tau1; tau2];
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 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
301
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
302 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
303 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
304
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
305 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
306
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
307 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
308 end
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
309
9a647dcccbdd Added pausing option to noname.animate. Added characteristic bc to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 54
diff changeset
310 closure = @closure_fun;
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
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
313
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 % 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
315 function closure = boundary_condition_inflow(obj, boundary)
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
316 [~,~,s] = obj.get_boundary_ops(boundary);
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
317
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
318 switch s
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 2];
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
321 case 1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
322 p_in = [1 3];
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
323 end
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
324
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
325 a = obj.gamma - 1;
63
21b18edb9667 Fixed error in bc.
Jonatan Werpers <jonatan@werpers.com>
parents: 62
diff changeset
326 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
327 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
328 0 -1/2*u*a a; %p
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
329 ];
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
330
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
331 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
332 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
333 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
334
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
335 % 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
336 function closure = boundary_condition_outflow(obj, boundary)
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
337 [~,~,s] = obj.get_boundary_ops(boundary);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
338
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
339 switch s
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
340 case -1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
341 p_in = 2;
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
342 case 1
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
343 p_in = 3;
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
344 end
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
345
60
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
346 a = obj.gamma -1;
e707cd9e6d0a Euler1d: Created a general bc routine.
Jonatan Werpers <jonatan@werpers.com>
parents: 59
diff changeset
347 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
348
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
349 closure = boundary_condition_L(obj, boundary, L, p_in);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
350 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
351
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
352 % 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
353 function closure = boundary_condition_inflow_rho(obj, boundary)
91
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
354 [~,~,s] = obj.get_boundary_ops(boundary);
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
355
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
356 switch s
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
357 case -1
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
358 p_in = [1 2];
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
359 case 1
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
360 p_in = [1 3];
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
361 end
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
362
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
363 a = obj.gamma - 1;
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
364 L = @(rho,~,~)[
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
365 0 1/rho 0;
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
366 1 0 0;
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
367 ];
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
368
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
369 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
370 end
2102af217134 Added inflow bc for rho
Jonatan Werpers <jonatan@werpers.com>
parents: 78
diff changeset
371
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
372 % 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
373 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
374 [~,~,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
375
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 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
377 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
378 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
379 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
380 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
381 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
382
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 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
384
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
385 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
386 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
387
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
388 % Set wall boundary condition v = 0.
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
389 function closure = boundary_condition_wall(obj,boundary)
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
390 [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
391
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
392 % 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
393 % vi låta den vara fri?
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
394
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 switch s
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
397 case -1
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
398 p_in = 2;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
399 p_zero = 1;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
400 p_ot = 3;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
401 case 1
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
402 p_in = 3;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
403 p_zero = 1;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
404 p_ot = 2;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
405 otherwise
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
406 error();
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
407 end
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
408
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
409 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
410 pt(p) = 1:length(p); % Inverse permutation
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
411
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
412 function o = closure_fun(q)
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
413
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
414 q_s = e_S'*q;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
415 rho = q_s(1);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
416 u = q_s(2)/rho;
50
75ebf5d3cfe5 Performance improvement for Euler1d
Jonatan Werpers <jonatan@werpers.com>
parents: 49
diff changeset
417 c = obj.c(q_s);
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
418
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
419 T = obj.T(q_s);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
420 R = -(u-c)/(u+c);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
421 % l = [u, u+c, u-c];
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 % p_in = find(s*l <= 0);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
424 % p_ot = find(s*l > 0);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
425
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 tau1 = -2*c;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
428 tau2 = [0; 0]; % Penalty only on ingoing char.
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
429
54
2194cd385419 Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
Jonatan Werpers <jonatan@werpers.com>
parents: 50
diff changeset
430 % 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
431 % Lambda_ot = diag(abs(l(p_ot)));
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
432
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
433 tauHat = [tau1; tau2];
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
434 tau = -s*e_S*sparse(T*tauHat(pt));
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 w_s = inv(T)*q_s;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
437 % w_s = T\q_s;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
438 % w_s = Tinv * q_s; % Med analytisk matris
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
439 w_in = w_s(p_in);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
440 w_ot = w_s(p_ot);
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
441
78
80948a4084f3 Added the factor or 1/2 to the wall bc term as well.
Jonatan Werpers <jonatan@werpers.com>
parents: 77
diff changeset
442 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
443 end
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 closure = @closure_fun;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
446 end
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
447
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
448 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
449 error('NOT DONE')
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
450 % u denotes the solution in the own domain
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
451 % v denotes the solution in the neighbour domain
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
452 [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
453 [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
454
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
455 tuning = 2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
456
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
457 alpha_u = obj.alpha;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
458 alpha_v = neighbour_scheme.alpha;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
459
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
460 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
461 % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
462 tau4 = s_u*alpha_u/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
463
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
464 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
465 sig3 = -s_u*alpha_u/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
466
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
467 phi2 = s_u*1/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
468
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
469 psi1 = -s_u*1/2;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
470
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
471 tau = tau1*e_u + tau4*d3_u;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
472 sig = sig2*d1_u + sig3*d2_u ;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
473 phi = phi2*d1_u ;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
474 psi = psi1*e_u ;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
475
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
476 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
477 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
478 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
479
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
480 % 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
481 % 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
482 function [e,E,s] = get_boundary_ops(obj,boundary)
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
483 switch boundary
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
484 case 'l'
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
485 e = obj.e_l;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
486 E = obj.e_L;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
487 s = -1;
49
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
488 case 'r'
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
489 e = obj.e_r;
8f0c2dc747dd Made lots of updates to Euler1d.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
490 E = obj.e_R;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
491 s = 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
492 otherwise
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
493 error('No such boundary: boundary = %s',boundary);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
494 end
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
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
497 function N = size(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
498 N = prod(obj.m);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
499 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
500
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
501 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
502 end