Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 49:8f0c2dc747dd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Made lots of updates to Euler1d.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 05 Nov 2015 16:33:53 -0800
parents 48b6fb693025
children 75ebf5d3cfe5
comparison
equal deleted inserted replaced
48:b21c53ff61d4 49:8f0c2dc747dd
1 classdef SchmBeam2d < noname.Scheme 1 classdef Euler1d < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 N % Number of points total 4 N % Number of points total
5 h % Grid spacing 5 h % Grid spacing
6 u % Grid values 6 u % Grid values
7 x % Values of x and y for each 7 x % Values of x and y for each
8 order % Order accuracy for the approximation 8 order % Order accuracy for the approximation
9 9
10 D % non-stabalized scheme operator 10 D % non-stabalized scheme operator
11 Fx
11 M % Derivative norm 12 M % Derivative norm
12 alpha 13 gamma
14
15 T
16 p
17
13 18
14 H % Discrete norm 19 H % Discrete norm
15 Hi 20 Hi
16 e_l, e_r 21 e_l, e_r, e_L, e_R;
17 22
18 end 23 end
19 24
20 methods 25 methods
21 function obj = SchmBeam2d(m,xlim,order,gamma,opsGen) 26 function obj = Euler1d(m,xlim,order,gama,opsGen,do_upwind)
22 default_arg('opsGen',@sbp.Ordinary); 27 default_arg('opsGen',@sbp.Ordinary);
23 default_arg('gamma', 1.4); 28 default_arg('gama', 1.4);
24 29 default_arg('do_upwind', false);
25 [x, h] = util.get_grid(xlim{:},m_x); 30 gamma = gama;
26 31
27 ops = opsGen(m_x,h_x,order); 32 [x, h] = util.get_grid(xlim{:},m);
28 33
29 I_x = speye(m); 34 if do_upwind
30 I_3 = speye(3); 35 ops = spb.Upwind(m,h,order);
31 36 Dp = ops.derivatives.Dp;
32 D1 = sparse(ops.derivatives.D1); 37 Dm = ops.derivatives.Dm;
38
39 printExpr('issparse(Dp)');
40 printExpr('issparse(Dm)');
41
42 D1 = (Dp + Dm)/2;
43 else
44 ops = opsGen(m,h,order);
45 D1 = sparse(ops.derivatives.D1);
46 end
47
33 H = sparse(ops.norms.H); 48 H = sparse(ops.norms.H);
34 Hi = sparse(ops.norms.HI); 49 Hi = sparse(ops.norms.HI);
35 e_l = sparse(ops.boundary.e_1); 50 e_l = sparse(ops.boundary.e_1);
36 e_r = sparse(ops.boundary.e_m); 51 e_r = sparse(ops.boundary.e_m);
37 52
53 I_x = speye(m);
54 I_3 = speye(3);
55
56
38 D1 = kr(D1, I_3); 57 D1 = kr(D1, I_3);
58 if do_upwind
59 Ddisp = kr(Ddisp,I_3);
60 end
39 61
40 % Norms 62 % Norms
41 obj.H = kr(H,I_3); 63 obj.H = kr(H,I_3);
64 obj.Hi = kr(Hi,I_3);
42 65
43 % Boundary operators 66 % Boundary operators
44 obj.e_l = kr(e_l,I_3); 67 obj.e_l = e_l;
45 obj.e_r = kr(e_r,I_3); 68 obj.e_r = e_r;
69 obj.e_L = kr(e_l,I_3);
70 obj.e_R = kr(e_r,I_3);
46 71
47 obj.m = m; 72 obj.m = m;
48 obj.h = h; 73 obj.h = h;
49 obj.order = order; 74 obj.order = order;
50 75
51
52 % Man har Q_t+F_x=0 i 1D Euler, där 76 % Man har Q_t+F_x=0 i 1D Euler, där
53 % q=[rho, rho*u, e]^T 77 % q=[rho, rho*u, e]^T
54 % F=[rho*u, rho*u^2+p, (e+p)*u] ^T 78 % F=[rho*u, rho*u^2+p, (e+p)*u] ^T
55 % p=(gamma-1)*(e-rho/2*u^2); 79 % p=(gamma-1)*(e-rho*u^2/2);
56
57 80
58 %Solving on form q_t + F_x = 0 81 %Solving on form q_t + F_x = 0
59 function o = F(q) 82 function o = F(q)
60 o = [q(2); q(2).^2/q(1) + p(q); (q(3)+p(q))*q(2)/q(1)]; 83 o = [q(2); q(2).^2/q(1) + p(q); (q(3)+p(q))*q(2)/q(1)];
61 end 84 end
62 85
63 % Equation of state 86 % Equation of state
64 function o = p(q) 87 function o = p(q)
65 o = (gamma-1)*(q(3)-q(2).^2/q(1)/2); 88 o = (gamma-1)*(q(3)-q(2).^2/q(1)/2);
66 end 89 end
67
68 90
69 % R = 91 % R =
70 % [sqrt(2*(gamma-1))*rho , rho , rho ; 92 % [sqrt(2*(gamma-1))*rho , rho , rho ;
71 % sqrt(2*(gamma-1))*rho*u , rho*(u+c) , rho*(u-c) ; 93 % sqrt(2*(gamma-1))*rho*u , rho*(u+c) , rho*(u-c) ;
72 % sqrt(2*(gamma-1))*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]); 94 % sqrt(2*(gamma-1))*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]);
73 function o = R(q) 95 function o = T(q)
74 rho = q(1); 96 rho = q(1);
75 u = q(2)/q(1); 97 u = q(2)/q(1);
76 e = q(3); 98 e = q(3);
99
100 c = sqrt(gamma*p(q)/rho);
77 101
78 sqrt2gamm = sqrt(2*(gamma-1)); 102 sqrt2gamm = sqrt(2*(gamma-1));
79 103
80 o = [ 104 o = [
81 sqrt2gamm*rho , rho , rho ; 105 sqrt2gamm*rho , rho , rho ;
82 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ; 106 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ;
83 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 107 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
84 ]; 108 ];
109 % Devide columns by stuff to get rid of extra comp?
85 end 110 end
86 111
87 function o = Fx(q) 112 function o = Fx(q)
88 o = zeros(size(q)); 113 o = zeros(size(q));
89 for i = 1:3:3*m 114 for i = 1:3:3*m
90 o(i:i+2) = F(q(i:i+2)); 115 o(i:i+2) = F(q(i:i+2));
91 end 116 end
92 end 117 o = D1*o;
93 118 end
94 119
95 120
96 % A=R*Lambda*inv(R), där Lambda=diag(u, u+c, u-c) (c är ljudhastigheten) 121 % A=R*Lambda*inv(R), där Lambda=diag(u, u+c, u-c) (c är ljudhastigheten)
97 % c^2=gamma*p/rho 122 % c^2=gamma*p/rho
98 % function o = A(rho,u,e) 123 % function o = A(rho,u,e)
99 % end 124 % end
100 125 if do_upwind
101 126 obj.D = @(q)-Fx(q) + Ddisp*(1)*u;
102 obj.D = @Fx; 127 else
128 obj.D = @(q)-Fx(q);
129 end
130
131 obj.Fx = @Fx;
132 obj.T = @T;
103 obj.u = x; 133 obj.u = x;
104 obj.x = kr(x,ones(3,1)); 134 obj.x = kr(x,ones(3,1));
105 end 135 obj.p = @p;
106 136 obj.gamma = gamma;
107 137 end
108 % Closure functions return the opertors applied to the own doamin to close the boundary 138
109 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 139
110 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 140 % Enforces the boundary conditions
111 % type is a string specifying the type of boundary condition if there are several. 141 % w+ = R*w- + g(t)
112 % data is a function returning the data that should be applied at the boundary. 142 function closure = boundary_condition(obj,boundary, type, varargin)
113 % neighbour_scheme is an instance of Scheme that should be interfaced to. 143 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
114 % neighbour_boundary is a string specifying which boundary to interface to.
115 function [closure, penalty] = boundary_condition(obj,boundary, alpha,data)
116 default_arg('alpha',0);
117 default_arg('data',0);
118 144
119 % Boundary condition on form 145 % Boundary condition on form
120 % w_in = w_out + g, where g is data 146 % w_in = R*w_out + g, where g is data
121 147
122 [e,s] = obj.get_boundary_ops(boundary); 148 % How to handle when the number of BC we want to set changes
123 149 % How to handel u = 0 as an example
124 tuning = 1; % ????????????????????????? 150
125 151 % Måste sätta in s och fixa tecken som de ska vara
126 tau = R(q)*lambda(q)*tuning; % SHOULD THIS BE abs(lambda)????? 152
127 153 % Kanske ska man tala om vilka karakteristikor som är in och ut i anropet
128 function closure_fun(q,t) 154 % Man kan sen kolla om det stämmer (men hur blir det med u=0?)
129 q_b = e * q; 155 % Och antingen tillåta att man skickar in flera alternativ och väljer automatiskt
130 end 156 % eller låta koden utanför bestämma vilken penalty som ska appliceras.
131 157
132 function penalty_fun(q,t) 158 switch type
133 end 159 case 'wall'
134 160 closure = obj.boundary_condition_wall(boundary);
135
136
137
138
139 % tau1 < -alpha^2/gamma
140
141 tau1 = tuning * alpha/delt;
142 tau4 = s*alpha;
143
144 sig2 = tuning * alpha/gamm;
145 sig3 = -s*alpha;
146
147 tau = tau1*e+tau4*d3;
148 sig = sig2*d1+sig3*d2;
149
150 closure = halfnorm_inv*(tau*e' + sig*d1');
151
152 pp_e = halfnorm_inv*tau;
153 pp_d = halfnorm_inv*sig;
154 switch class(data)
155 case 'double'
156 penalty_e = pp_e*data;
157 penalty_d = pp_d*data;
158 case 'function_handle'
159 penalty_e = @(t)pp_e*data(t);
160 penalty_d = @(t)pp_d*data(t);
161 otherwise 161 otherwise
162 error('Wierd data argument!') 162 error('Unsupported bc type: %s', type);
163 end 163 end
164 164
165 % T = obj.T(e_S*v);
166
167 % c = sqrt(gamma*p/rho);
168 % l = [u, u+c, u-c];
169
170 % p_in = find(s*l <= 0);
171 % p_ot = find(s*l > 0);
172
173 % p = [p_in, p_ot] % Permutation to sort
174 % pt(p) = 1:length(p); % Inverse permutation
175
176 % L_in = diag(abs(l(p_in)));
177 % L_ot = diag(abs(l(p_ot)));
178
179 % Lambda = diag(u, u+c, u-c);
180
181 % tau = -e_S*sparse(T*[L_in; R'*L_in]); % Något med pt
182
183
184 % w_s = T'*(e_S*v);
185 % w_in = w_s(p_in);
186 % w_ot = w_s(p_ot);
187
188
189 % function closure_fun(q,t)
190 % obj.Hi * tau * (w_in - R*w_ot - g(t));
191 % end
192
193 % function closure_fun_indep(q,t)
194 % obj.Hi * tau * (w_in - R*w_ot - g;
195 % end
196
197
198 % switch class(g)
199 % case 'double'
200 % closure = @closure_fun;
201 % case 'function_handle'
202 % closure = @closure_fun_indep;
203 % otherwise
204 % error('Wierd data argument!');
205 % end
206
207 end
208
209 % Set wall boundary condition v = 0.
210 function closure = boundary_condition_wall(obj,boundary)
211 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
212
213 % v = 0 corresponds to
214 % L = [0 1 0];
215 % g = 0
216 %
217 % Tp =
218 % R = -(u-c)/(u+c)
219
220 % tau = alpha * (u+c)
221 % (alpha+1)(u+c) + 1/4* alpha^2|u-c| <= 0
222 % 4*(alpha+1)(u+c) + alpha^2|u-c| <= 0
223 % 4 * (alpha+1)(u+c) + alpha^2|u| + alpha^2*c <= 0
224 % |u|*(sgn(u)*4 + sgn(u)*4*alpha + alpha^2) + c*(4 + 4alpha + alpha^2) <= 0
225 % |u|*(alpha^2 + 4*sgn(u)*alpha + 4*sgn(u)) + c*(alpha+2)^2 <= 0
226 % |u|*[(alpha + 2*sgn(u))^2 - 4*(sgn(u)-1)] + c*(alpha+2)^2 <= 0
227
228
229 % om vi låtsas att u = 0:
230 % (alpha+1)c + 1/4 * alpha^2*c <= 0
231 % alpha^2 + 4*alpha +4 <= 0
232 % (alpha + 2)^2 <= 0
233 % alpha = -2 gives tau = -2*c;
234
235
236 % Vill vi sätta penalty på karateristikan som är nära noll också eller vill
237 % vi låta den vara fri?
238
239
240 switch s
241 case -1
242 p_in = 2;
243 p_zero = 1;
244 p_ot = 3;
245 case 1
246 p_in = 3;
247 p_zero = 1;
248 p_ot = 2;
249 otherwise
250 error();
251 end
252
253 p = [p_in, p_zero, p_ot]; % Permutation to sort
254 pt(p) = 1:length(p); % Inverse permutation
255
256 function o = closure_fun(q)
257 p = obj.p(q);
258
259 q_s = e_S'*q;
260 rho = q_s(1);
261 u = q_s(2)/rho;
262 c = sqrt(obj.gamma*p/rho);
263
264 T = obj.T(q_s);
265 R = -(u-c)/(u+c);
266 % l = [u, u+c, u-c];
267
268 % p_in = find(s*l <= 0);
269 % p_ot = find(s*l > 0);
270
271
272 tau1 = -2*c;
273 tau2 = [0; 0]; % Penalty only on ingoing char.
274
275 % L_in = diag(abs(l(p_in)));
276 % L_ot = diag(abs(l(p_ot)));
277
278 tauHat = [tau1; tau2];
279 tau = -s*e_S*sparse(T*tauHat(pt));
280
281 w_s = inv(T)*q_s;
282 % w_s = T\q_s;
283 % w_s = Tinv * q_s; % Med analytisk matris
284 w_in = w_s(p_in);
285 w_ot = w_s(p_ot);
286
287 o = obj.Hi * tau * (w_in - R*w_ot);
288 end
289
290 closure = @closure_fun;
165 end 291 end
166 292
167 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 293 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
294 error('NOT DONE')
168 % u denotes the solution in the own domain 295 % u denotes the solution in the own domain
169 % v denotes the solution in the neighbour domain 296 % v denotes the solution in the neighbour domain
170 [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); 297 [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
171 [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 298 [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
172 299
195 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); 322 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
196 end 323 end
197 324
198 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 325 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
199 % The right boundary is considered the positive boundary 326 % The right boundary is considered the positive boundary
200 function [e,d1,d2,d3,s,gamm, delt, halfnorm_inv] = get_boundary_ops(obj,boundary) 327 function [e,E,s] = get_boundary_ops(obj,boundary)
201 switch boundary 328 switch boundary
202 case 'w' 329 case 'l'
203 e = obj.e_w; 330 e = obj.e_l;
204 d1 = obj.d1_w; 331 E = obj.e_L;
205 d2 = obj.d2_w;
206 d3 = obj.d3_w;
207 s = -1; 332 s = -1;
208 gamm = obj.gamm_x; 333 case 'r'
209 delt = obj.delt_x; 334 e = obj.e_r;
210 halfnorm_inv = obj.Hix; 335 E = obj.e_R;
211 case 'e'
212 e = obj.e_e;
213 d1 = obj.d1_e;
214 d2 = obj.d2_e;
215 d3 = obj.d3_e;
216 s = 1; 336 s = 1;
217 gamm = obj.gamm_x;
218 delt = obj.delt_x;
219 halfnorm_inv = obj.Hix;
220 case 's'
221 e = obj.e_s;
222 d1 = obj.d1_s;
223 d2 = obj.d2_s;
224 d3 = obj.d3_s;
225 s = -1;
226 gamm = obj.gamm_y;
227 delt = obj.delt_y;
228 halfnorm_inv = obj.Hiy;
229 case 'n'
230 e = obj.e_n;
231 d1 = obj.d1_n;
232 d2 = obj.d2_n;
233 d3 = obj.d3_n;
234 s = 1;
235 gamm = obj.gamm_y;
236 delt = obj.delt_y;
237 halfnorm_inv = obj.Hiy;
238 otherwise 337 otherwise
239 error('No such boundary: boundary = %s',boundary); 338 error('No such boundary: boundary = %s',boundary);
240 end 339 end
241 end 340 end
242 341