Mercurial > repos > public > sbplib

annotate +scheme/Euler1d.m @ 59:e431c1260f52

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