Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 77:0b07ff8a0a12

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed bugs and added outflow bc for rho. Added missing abs() added missing 1/2.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 25 Nov 2015 15:19:55 +0100
parents 21b18edb9667
children 80948a4084f3
comparison
equal deleted inserted replaced
76:5c569cbef49e 77:0b07ff8a0a12
168 closure = obj.boundary_condition_wall(boundary); 168 closure = obj.boundary_condition_wall(boundary);
169 case 'inflow' 169 case 'inflow'
170 closure = obj.boundary_condition_inflow(boundary,varargin{:}); 170 closure = obj.boundary_condition_inflow(boundary,varargin{:});
171 case 'outflow' 171 case 'outflow'
172 closure = obj.boundary_condition_outflow(boundary,varargin{:}); 172 closure = obj.boundary_condition_outflow(boundary,varargin{:});
173 case 'outflow_rho'
174 closure = obj.boundary_condition_outflow_rho(boundary,varargin{:});
173 case 'char' 175 case 'char'
174 closure = obj.boundary_condition_char(boundary,varargin{:}); 176 closure = obj.boundary_condition_char(boundary,varargin{:});
175 otherwise 177 otherwise
176 error('Unsupported bc type: %s', type); 178 error('Unsupported bc type: %s', type);
177 end 179 end
205 Tin = T(:,p_in); 207 Tin = T(:,p_in);
206 Tot = T(:,p_ot); 208 Tot = T(:,p_ot);
207 209
208 % Convert bc from ordinary form to characteristic form. 210 % Convert bc from ordinary form to characteristic form.
209 % Lq = g => w_in = Rw_ot + g_tilde 211 % Lq = g => w_in = Rw_ot + g_tilde
212
213 Lambda = obj.Lambda(q_s);
214
215 % Setup the penalty parameter
216 tau1 = -2*abs(Lambda(p_in,p_in));
217 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
218
219 tauHat = [tau1; tau2];
220 tau = -s*e_S*sparse(T*tauHat(pt,:));
221
210 L = L_fun(rho,u,e); 222 L = L_fun(rho,u,e);
211 g = g_fun(t); 223 g = g_fun(t);
212 224
213 Lambda = obj.Lambda(q_s); 225 % printExpr('s')
214 226 % penalty = tauHat(pt,:)*inv(L*Tin)*(L*q_s - g);
215 R =-inv(L*Tin)*(L*Tot); 227 % tauHatPt = tauHat(pt,:);
216 g_tilde = inv(L*Tin)*g; 228 % display(tauHatPt);
217 229 % display(penalty);
218 % Setup the penalty parameter 230 % pause
219 tau1 = -2*Lambda(p_in,p_in); 231
220 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char. 232 o = 1/2*obj.Hi * tau * inv(L*Tin)*(L*q_s - g);
221
222 tauHat = [tau1; tau2];
223 tau = -s*e_S*sparse(T*tauHat(pt,:));
224
225 % Calculate the numerical caracteristics and apply the penalty
226 w_s = inv(T)*q_s;
227 % w_s = T\q_s;
228 % w_s = Tinv * q_s; % Med analytisk matris
229 w_in = w_s(p_in);
230 w_ot = w_s(p_ot);
231
232 o = obj.Hi * tau * (w_in - R*w_ot - g_tilde);
233 end 233 end
234 closure = @closure_fun; 234 closure = @closure_fun;
235 end 235 end
236 236
237 function closure = boundary_condition_char(obj,boundary,w_data) 237 function closure = boundary_condition_char(obj,boundary,w_data)
252 p = [p_in p_ot]; 252 p = [p_in p_ot];
253 pt(p) = 1:length(p); 253 pt(p) = 1:length(p);
254 254
255 T = obj.T(q_s); 255 T = obj.T(q_s);
256 256
257 tau1 = -2*diag(Lambda(p_in)); 257 tau1 = -2*diag(abs(Lambda(p_in)));
258 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char. 258 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
259 259
260 tauHat = [tau1; tau2]; 260 tauHat = [tau1; tau2];
261 261
262 tau = -s*e_S*sparse(T*tauHat(pt,:)); 262 tau = -s*e_S*sparse(T*tauHat(pt,:));
265 w_in = w_s(p_in); 265 w_in = w_s(p_in);
266 266
267 w_s_data = w_data(t); 267 w_s_data = w_data(t);
268 w_in_data = w_s_data(p_in); 268 w_in_data = w_s_data(p_in);
269 269
270 o = obj.Hi * tau * (w_in - w_in_data); 270 o = 1/2*obj.Hi * tau * (w_in - w_in_data);
271 end 271 end
272 272
273 closure = @closure_fun; 273 closure = @closure_fun;
274 274
275 end 275 end
309 end 309 end
310 310
311 a = obj.gamma -1; 311 a = obj.gamma -1;
312 L = @(~,u,~)a*[0 -1/2*u 1]; 312 L = @(~,u,~)a*[0 -1/2*u 1];
313 g = @(t)[p_data(t)]; 313 g = @(t)[p_data(t)];
314
315
316 closure = boundary_condition_L(obj, boundary, L, g, p_in);
317
318 end
319
320 function closure = boundary_condition_outflow_rho(obj, boundary, rho_data)
321 [~,~,s] = obj.get_boundary_ops(boundary);
322
323 switch s
324 case -1
325 p_in = 2;
326 case 1
327 p_in = 3;
328 end
329
330 L = @(~,~,~)[1 0 0];
331 g = @(t)[rho_data(t)];
314 332
315 333
316 closure = boundary_condition_L(obj, boundary, L, g, p_in); 334 closure = boundary_condition_L(obj, boundary, L, g, p_in);
317 335
318 end 336 end