Mercurial > repos > public > sbplib
changeset 54:2194cd385419
Euler1d: Fixed russian dissipation. Cleaned up. Added non-working BC.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 11 Nov 2015 19:21:38 -0800 |
| parents | cca09cc5121d |
| children | a8ed986fcf57 |
| files | +scheme/Euler1d.m |
| diffstat | 1 files changed, 137 insertions(+), 59 deletions(-) [+] |
line wrap: on
line diff
diff -r cca09cc5121d -r 2194cd385419 +scheme/Euler1d.m --- a/+scheme/Euler1d.m Wed Nov 11 19:20:32 2015 -0800 +++ b/+scheme/Euler1d.m Wed Nov 11 19:21:38 2015 -0800 @@ -33,14 +33,12 @@ [x, h] = util.get_grid(xlim{:},m); if do_upwind - ops = spb.Upwind(m,h,order); + ops = sbp.Upwind(m,h,order); Dp = ops.derivatives.Dp; Dm = ops.derivatives.Dm; - printExpr('issparse(Dp)'); - printExpr('issparse(Dm)'); - D1 = (Dp + Dm)/2; + Ddisp = (Dp - Dm)/2; else ops = opsGen(m,h,order); D1 = sparse(ops.derivatives.D1); @@ -54,7 +52,6 @@ I_x = speye(m); I_3 = speye(3); - D1 = kr(D1, I_3); if do_upwind Ddisp = kr(Ddisp,I_3); @@ -96,12 +93,9 @@ o = sqrt(gamma*p(Q)./Q(1,:)); end - - % R = - % [sqrt(2*(gamma-1))*rho , rho , rho ; - % sqrt(2*(gamma-1))*rho*u , rho*(u+c) , rho*(u-c) ; - % 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]); function o = T(q) + % T is the transformation such that A = T*Lambda*inv(T) + % where Lambda=diag(u, u+c, u-c) rho = q(1); u = q(2)/q(1); e = q(3); @@ -113,9 +107,9 @@ o = [ sqrt2gamm*rho , rho , rho ; sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ; - 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 - ]; - % Devide columns by stuff to get rid of extra comp? + 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 ; + ]; + % Devide columns by stuff to get rid of extra comp? end function o = Fx(q) @@ -124,13 +118,25 @@ o = D1*o; end + function o = Fx_disp(q) + Q = reshape(q,3,m); + f = reshape(F(Q),3*m,1); - % A=R*Lambda*inv(R), där Lambda=diag(u, u+c, u-c) (c är ljudhastigheten) - % c^2=gamma*p/rho - % function o = A(rho,u,e) - % end + c = c(Q); + lambda_max = c+abs(Q(2,:)./Q(1,:)); + % lambda_max = max(lambda_max); + + lamb_Q(1,:) = lambda_max.*Q(1,:); + lamb_Q(2,:) = lambda_max.*Q(2,:); + lamb_Q(3,:) = lambda_max.*Q(3,:); + + lamb_q = reshape(lamb_Q,3*m, 1); + + o = -D1*f + Ddisp*lamb_q; + end + if do_upwind - obj.D = @(q)-Fx(q) + Ddisp*(1)*u; + obj.D = @Fx_disp; else obj.D = @(q)-Fx(q); end @@ -153,64 +159,136 @@ % Boundary condition on form % w_in = R*w_out + g, where g is data - % How to handle when the number of BC we want to set changes - % How to handel u = 0 as an example - - % Måste sätta in s och fixa tecken som de ska vara - - % Kanske ska man tala om vilka karakteristikor som är in och ut i anropet - % Man kan sen kolla om det stämmer (men hur blir det med u=0?) - % Och antingen tillåta att man skickar in flera alternativ och väljer automatiskt - % eller låta koden utanför bestämma vilken penalty som ska appliceras. - switch type case 'wall' closure = obj.boundary_condition_wall(boundary); + case 'inflow' + closure = obj.boundary_condition_inflow(boundary,varargin{:}); + case 'outflow' + closure = obj.boundary_condition_outflow(boundary,varargin{:}); otherwise error('Unsupported bc type: %s', type); end - % T = obj.T(e_S*v); + end + + function closure = boundary_condition_inflow(obj, boundary, p_data, v_data) + [e_s,e_S,s] = obj.get_boundary_ops(boundary); - % c = sqrt(gamma*p/rho); - % l = [u, u+c, u-c]; + switch s + case -1 + p_in = [1 2]; + p_ot = [3]; + case 1 + p_in = [1 3]; + p_ot = [2]; + otherwise + error(); + end - % p_in = find(s*l <= 0); - % p_ot = find(s*l > 0); + p = [p_in, p_ot]; % Permutation to sort + pt(p) = 1:length(p); % Inverse permutation + + gamma = obj.gamma; - % p = [p_in, p_ot] % Permutation to sort - % pt(p) = 1:length(p); % Inverse permutation + function o = closure_fun(q,t) + q_s = e_S'*q; + rho = q_s(1); + u = q_s(2)/rho; + e = q_s(3); + + c = obj.c(q_s); + + T = obj.T(q_s); - % L_in = diag(abs(l(p_in))); - % L_ot = diag(abs(l(p_ot))); + % L = [ + % 0 1 0; + % 0 -1/2*u 1; + % ]; + LTpi = [ + 2^(1/2)/(2*rho*u*(gamma - 1)^(1/2)), -(2^(1/2)*(u - c*s))/(u*(2*e*gamma - rho*(gamma*u^2 - c*s*u + 2*c*s))*(gamma - 1)^(1/2)); + 0, 2/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2)); + ]; + R = [ + (2^(1/2)*c*s*(rho*(gamma - 1)*u^2 + 2*rho*u - 2*e*gamma))/(u*(gamma - 1)^(1/2)*(- gamma*rho*u^2 + c*rho*s*u + 2*e*gamma - 2*c*rho*s)); + (2*c*rho*s*(u - 2))/(2*e*gamma - gamma*rho*u^2 + c*rho*s*(u - 2)) - 1; + ]; - % Lambda = diag(u, u+c, u-c); - - % tau = -e_S*sparse(T*[L_in; R'*L_in]); % Något med pt + tau1 = -2*[ + u, 0; + 0, c + abs(u); + ]; + tau2 = [0 0]; % Penalty only on ingoing char. - % w_s = T'*(e_S*v); - % w_in = w_s(p_in); - % w_ot = w_s(p_ot); + tauHat = [tau1; tau2]; + tau = -s*e_S*sparse(T*tauHat(pt,:)); + + w_s = inv(T)*q_s; + % w_s = T\q_s; + % w_s = Tinv * q_s; % Med analytisk matris + w_in = w_s(p_in); + w_ot = w_s(p_ot); + + + o = obj.Hi * tau * (w_in - R*w_ot - LTpi*[v_data(t); p_data(t)]); + end + + closure = @closure_fun; + end - % function closure_fun(q,t) - % obj.Hi * tau * (w_in - R*w_ot - g(t)); - % end + function closure = boundary_condition_outflow(obj, boundary, p_data) + [e_s,e_S,s] = obj.get_boundary_ops(boundary); - % function closure_fun_indep(q,t) - % obj.Hi * tau * (w_in - R*w_ot - g; - % end + switch s + case -1 + p_in = 2; + p_ot = [1 3]; + case 1 + p_in = 3; + p_ot = [1 2]; + otherwise + error(); + end + + p = [p_in, p_ot]; % Permutation to sort + pt(p) = 1:length(p); % Inverse permutation + + gamma = obj.gamma; + + function o = closure_fun(q,t) + q_s = e_S'*q; + rho = q_s(1); + u = q_s(2)/rho; + e = q_s(3); + + c = obj.c(q_s); - % switch class(g) - % case 'double' - % closure = @closure_fun; - % case 'function_handle' - % closure = @closure_fun_indep; - % otherwise - % error('Wierd data argument!'); - % end + T = obj.T(q_s); + + L = (gamma -1)*[0 -1/2*u 1]; + LTpi = 1/(L*T(:,p_in)); + R = -[0 (gamma*(e - rho*u^2/2) - rho*u*c)/(gamma*(e - rho*u^2/2) + rho*u*c)]; + + tau1 = [-2*(c - abs(u))]; + tau2 = [0; 0]; % Penalty only on ingoing char. + + + tauHat = [tau1; tau2]; + tau = -s*e_S*sparse(T*tauHat(pt)); + + w_s = inv(T)*q_s; + % w_s = T\q_s; + % w_s = Tinv * q_s; % Med analytisk matris + w_in = w_s(p_in); + w_ot = w_s(p_ot); + + o = obj.Hi * tau * (w_in - R*w_ot - LTpi*p_data(t)); + end + + closure = @closure_fun; end @@ -279,8 +357,8 @@ tau1 = -2*c; tau2 = [0; 0]; % Penalty only on ingoing char. - % L_in = diag(abs(l(p_in))); - % L_ot = diag(abs(l(p_ot))); + % Lambda_in = diag(abs(l(p_in))); + % Lambda_ot = diag(abs(l(p_ot))); tauHat = [tau1; tau2]; tau = -s*e_S*sparse(T*tauHat(pt));