Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 57:9a647dcccbdd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added pausing option to noname.animate. Added characteristic bc to Euler1d.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 12 Nov 2015 17:31:33 -0800
parents 2194cd385419
children 24103284e09d
comparison
equal deleted inserted replaced
56:ce90abc350c5 57:9a647dcccbdd
164 closure = obj.boundary_condition_wall(boundary); 164 closure = obj.boundary_condition_wall(boundary);
165 case 'inflow' 165 case 'inflow'
166 closure = obj.boundary_condition_inflow(boundary,varargin{:}); 166 closure = obj.boundary_condition_inflow(boundary,varargin{:});
167 case 'outflow' 167 case 'outflow'
168 closure = obj.boundary_condition_outflow(boundary,varargin{:}); 168 closure = obj.boundary_condition_outflow(boundary,varargin{:});
169 case 'char'
170 closure = obj.boundary_condition_char(boundary,varargin{:});
169 otherwise 171 otherwise
170 error('Unsupported bc type: %s', type); 172 error('Unsupported bc type: %s', type);
171 end 173 end
174
175 end
176
177 function closure = boundary_condition_char(obj,boundary,w_data)
178 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
179
180 function o = closure_fun(q,t)
181 q_s = e_S'*q;
182 rho = q_s(1);
183 u = q_s(2)/rho;
184 e = q_s(3);
185
186 c = obj.c(q_s);
187
188 Lambda = [u, u+c, u-c];
189
190 p_in = find(s*Lambda < 0);
191 p_ot = find(s*Lambda >= 0);
192 p = [p_in p_ot];
193 pt(p) = 1:length(p);
194
195 T = obj.T(q_s);
196
197 tau1 = -2*diag(Lambda(p_in));
198 tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
199
200 tauHat = [tau1; tau2];
201
202 tau = -s*e_S*sparse(T*tauHat(pt,:));
203
204 w_s = inv(T)*q_s;
205 w_in = w_s(p_in);
206
207 w_s_data = w_data(t);
208 w_in_data = w_s_data(p_in);
209
210 o = obj.Hi * tau * (w_in - w_in_data);
211 end
212
213 closure = @closure_fun;
172 214
173 end 215 end
174 216
175 function closure = boundary_condition_inflow(obj, boundary, p_data, v_data) 217 function closure = boundary_condition_inflow(obj, boundary, p_data, v_data)
176 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 218 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
220 ]; 262 ];
221 tau2 = [0 0]; % Penalty only on ingoing char. 263 tau2 = [0 0]; % Penalty only on ingoing char.
222 264
223 265
224 tauHat = [tau1; tau2]; 266 tauHat = [tau1; tau2];
267
225 tau = -s*e_S*sparse(T*tauHat(pt,:)); 268 tau = -s*e_S*sparse(T*tauHat(pt,:));
226 269
227 w_s = inv(T)*q_s; 270 w_s = inv(T)*q_s;
228 % w_s = T\q_s; 271 % w_s = T\q_s;
229 % w_s = Tinv * q_s; % Med analytisk matris 272 % w_s = Tinv * q_s; % Med analytisk matris