Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 97:33dba20b4b9d feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged with deafult.
author Martin Almquist <martin.almquist@it.uu.se>
date Wed, 02 Dec 2015 16:29:54 +0100
parents 5c41941b9e5e
children ce4eecbcb915
comparison
equal deleted inserted replaced
88:7b092f0fd620 97:33dba20b4b9d
15 Hi 15 Hi
16 e_l, e_r, e_L, e_R; 16 e_l, e_r, e_L, e_R;
17 17
18 end 18 end
19 19
20 properties (Constant)
21 SUBSONIC_INFLOW = 1;
22 SUBSONIC_OUTFLOW = -1;
23 NO_FLOW = 0;
24 SUPERSONIC_INFLOW = 2;
25 SUPERSONIC_OUTFLOW = -2;
26 end
27
20 methods 28 methods
21 function obj = Euler1d(m,xlim,order,gama,opsGen,do_upwind) 29 function obj = Euler1d(m,xlim,order,gama,opsGen,do_upwind)
22 default_arg('opsGen',@sbp.Ordinary); 30 default_arg('opsGen',@sbp.Ordinary);
23 default_arg('gama', 1.4); 31 default_arg('gama', 1.4);
24 default_arg('do_upwind', false); 32 default_arg('do_upwind', false);
153 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 ; 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 ;
154 ]; 162 ];
155 % Devide columns by stuff to get rid of extra comp? 163 % Devide columns by stuff to get rid of extra comp?
156 end 164 end
157 165
166 function fs = flowStateL(obj, q)
167 q_l = obj.e_L'*q;
168 c = obj.c(q);
169 v = q_l(2,:)/q_l(1,:);
170
171 if v > c
172 fs = scheme.Euler1d.SUPERSONIC_INFLOW;
173 elseif v > 0
174 fs = scheme.Euler1d.SUBSONIC_INFLOW;
175 elseif v > -c
176 fs = scheme.Euler1d.SUBSONIC_OUTFLOW;
177 else
178 fs = scheme.Euler1d.SUPERSONIC_INFLOW;
179 end
180 end
181
182 % returns positiv values for inlfow, negative for outflow.
183 % +-1 for subsonic
184 function fs = flowStateR(obj, q)
185 q_r = obj.e_R'*q;
186 c = obj.c(q);
187 v = q_r(2,:)/q_r(1,:);
188
189 if v < -c
190 fs = scheme.Euler1d.SUPERSONIC_INFLOW;
191 elseif v < 0
192 fs = scheme.Euler1d.SUBSONIC_INFLOW;
193 elseif v < c
194 fs = scheme.Euler1d.SUBSONIC_OUTFLOW;
195 else
196 fs = scheme.Euler1d.SUPERSONIC_INFLOW;
197 end
198 end
199
158 % Enforces the boundary conditions 200 % Enforces the boundary conditions
159 % w+ = R*w- + g(t) 201 % w+ = R*w- + g(t)
160 function closure = boundary_condition(obj,boundary, type, varargin) 202 function closure = boundary_condition(obj,boundary, type, varargin)
161 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 203 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
162 204
168 closure = obj.boundary_condition_wall(boundary); 210 closure = obj.boundary_condition_wall(boundary);
169 case 'inflow' 211 case 'inflow'
170 closure = obj.boundary_condition_inflow(boundary,varargin{:}); 212 closure = obj.boundary_condition_inflow(boundary,varargin{:});
171 case 'outflow' 213 case 'outflow'
172 closure = obj.boundary_condition_outflow(boundary,varargin{:}); 214 closure = obj.boundary_condition_outflow(boundary,varargin{:});
173 case 'outflow_rho' 215 case 'inflow_rho'
216 closure = obj.boundary_condition_inflow_rho(boundary,varargin{:});
217 case 'outflow_rho'
174 closure = obj.boundary_condition_outflow_rho(boundary,varargin{:}); 218 closure = obj.boundary_condition_outflow_rho(boundary,varargin{:});
175 case 'char' 219 case 'char'
176 closure = obj.boundary_condition_char(boundary,varargin{:}); 220 closure = obj.boundary_condition_char(boundary,varargin{:});
177 otherwise 221 otherwise
178 error('Unsupported bc type: %s', type); 222 error('Unsupported bc type: %s', type);
205 % Calculate transformation matrix 249 % Calculate transformation matrix
206 T = obj.T(q_s); 250 T = obj.T(q_s);
207 Tin = T(:,p_in); 251 Tin = T(:,p_in);
208 Tot = T(:,p_ot); 252 Tot = T(:,p_ot);
209 253
210 % Convert bc from ordinary form to characteristic form. 254 % Calculate eigen value matrix
211 % Lq = g => w_in = Rw_ot + g_tilde
212
213 Lambda = obj.Lambda(q_s); 255 Lambda = obj.Lambda(q_s);
214 256
215 % Setup the penalty parameter 257 % Setup the penalty parameter
216 tau1 = -2*abs(Lambda(p_in,p_in)); 258 tau1 = -2*abs(Lambda(p_in,p_in));
217 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char. 259 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
218 260
219 tauHat = [tau1; tau2]; 261 tauHat = [tau1; tau2];
220 tau = -s*e_S*sparse(T*tauHat(pt,:)); 262 tau = e_S*sparse(T*tauHat(pt,:));
221 263
222 L = L_fun(rho,u,e); 264 L = L_fun(rho,u,e);
223 g = g_fun(t); 265 g = g_fun(t);
224
225 % printExpr('s')
226 % penalty = tauHat(pt,:)*inv(L*Tin)*(L*q_s - g);
227 % tauHatPt = tauHat(pt,:);
228 % display(tauHatPt);
229 % display(penalty);
230 % pause
231 266
232 o = 1/2*obj.Hi * tau * inv(L*Tin)*(L*q_s - g); 267 o = 1/2*obj.Hi * tau * inv(L*Tin)*(L*q_s - g);
233 end 268 end
234 closure = @closure_fun; 269 closure = @closure_fun;
235 end 270 end
313 g = @(t)[p_data(t)]; 348 g = @(t)[p_data(t)];
314 349
315 350
316 closure = boundary_condition_L(obj, boundary, L, g, p_in); 351 closure = boundary_condition_L(obj, boundary, L, g, p_in);
317 352
353 end
354
355 function closure = boundary_condition_inflow_rho(obj, boundary, rho_data, v_data)
356 [~,~,s] = obj.get_boundary_ops(boundary);
357
358 switch s
359 case -1
360 p_in = [1 2];
361 case 1
362 p_in = [1 3];
363 end
364
365 a = obj.gamma - 1;
366 L = @(rho,~,~)[
367 0 1/rho 0;
368 1 0 0;
369 ];
370 g = @(t)[
371 v_data(t);
372 rho_data(t);
373 ];
374
375 closure = boundary_condition_L(obj, boundary, L, g, p_in);
318 end 376 end
319 377
320 function closure = boundary_condition_outflow_rho(obj, boundary, rho_data) 378 function closure = boundary_condition_outflow_rho(obj, boundary, rho_data)
321 [~,~,s] = obj.get_boundary_ops(boundary); 379 [~,~,s] = obj.get_boundary_ops(boundary);
322 380