sbplib: +scheme/Euler1d.m comparison

comparison +scheme/Euler1d.m @ 60:e707cd9e6d0a

Euler1d: Created a general bc routine.

author	Jonatan Werpers <jonatan@werpers.com>
date	Sat, 14 Nov 2015 15:43:01 -0800
parents	e431c1260f52
children	183d9349b4c1

comparison

equal deleted inserted replaced

-:e431c1260f52
+:e707cd9e6d0a
 u % Grid values
 x % Values of x and y for each
 order % Order accuracy for the approximation
 D % non-stabalized scheme operator
-Fx
 M % Derivative norm
 gamma
 T
 p
 c
+Lambda
 H % Discrete norm
 Hi
 e_l, e_r, e_L, e_R;
 % Equation of state
 function o = p(Q)
 % Pressure p = (gamma-1)*(q3-q2.^2/q1/2)
 o = (gamma-1)*( Q(3,:)-1/2*Q(2,:).^2./Q(1,:) );
 end
+obj.p = @p;
 function o = c(Q)
 % Speed of light c = sqrt(obj.gamma*p/rho);
 o = sqrt(gamma*p(Q)./Q(1,:));
 end
+obj.c = @c;
+function o = Lambda(q)
+u = q(2)/q(1);
+c = obj.c(q);
+L = [u, u+c, u-c];
+o = diag(L);
+end
+obj.Lambda = @Lambda;
 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);
 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?
 end
+obj.T = @T;
 function o = Fx(q)
 Q = reshape(q,3,m);
 o = reshape(F(Q),3*m,1);
 o = D1*o;
 obj.D = @Fx_disp;
 else
 obj.D = @(q)-Fx(q);
 end
-obj.Fx = @Fx;
-obj.T = @T;
 obj.u = x;
 obj.x = kr(x,ones(3,1));
-obj.p = @p;
-obj.c = @c;
 obj.gamma = gamma;
 end
 % Enforces the boundary conditions
 error('Unsupported bc type: %s', type);
 end
 end
-function closure = boundary_condition_char(obj,boundary,w_data)
+% Sets the boundary condition Lq = g, where
+%   L = L(rho, u, e), g = g(t)
+%   p_in are the indecies of the ingoing characteristics.
+function closure = boundary_condition_L(obj, boundary, L_fun, g_fun, p_in)
 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
-function o = closure_fun(q,t)
+p_ot = 1:3;
-q_s = e_S'*q;
+p_ot(p_in) = [];
-rho = q_s(1);
-u = q_s(2)/rho;
-e = q_s(3);
-c = obj.c(q_s);
-Lambda = [u, u+c, u-c];
-p_in = find(s*Lambda < 0);
-p_ot = find(s*Lambda >= 0);
-p = [p_in p_ot];
-pt(p) = 1:length(p);
-T = obj.T(q_s);
-tau1 = -2*diag(Lambda(p_in));
-tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
-tauHat = [tau1; tau2];
-tau = -s*e_S*sparse(T*tauHat(pt,:));
-w_s = inv(T)*q_s;
-w_in = w_s(p_in);
-w_s_data = w_data(t);
-w_in_data = w_s_data(p_in);
-o = obj.Hi * tau * (w_in - w_in_data);
-end
-closure = @closure_fun;
-end
-function closure = boundary_condition_inflow(obj, boundary, p_data, v_data)
-[e_s,e_S,s] = obj.get_boundary_ops(boundary);
-switch s
-case -1
-p_in = [1 2];
-p_ot = [3];
-case 1
-p_in = [1 3];
-p_ot = [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)
 % Extract solution at the boundary
 q_s = e_S'*q;
 rho = q_s(1);
 Tin = T(:,p_in);
 Tot = T(:,p_ot);
 % Convert bc from ordinary form to characteristic form.
 %   Lq = g  => w_in = Rw_ot + g_tilde
-L = [
+L = L_fun(rho,u,e);
-0      1 0;
+g = g_fun(t);
-0 -1/2*u 1;
-];
+Lambda = obj.Lambda(q_s);
-g = [
-v_data(t);
-p_data(t);
-];
 R =-inv(L*Tin)*(L*Tot);
 g_tilde = inv(L*Tin)*g;
 % Setup the penalty parameter
-tau1 = -2*[
+tau1 = -2*Lambda(p_in,p_in);
-u,          0;
+tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
-0, c + abs(u);
-];
-tau2 = [0 0]; % Penalty only on ingoing char.
 tauHat = [tau1; tau2];
 tau = -s*e_S*sparse(T*tauHat(pt,:));
 % Calculate the numerical caracteristics and apply the penalty
 w_in = w_s(p_in);
 w_ot = w_s(p_ot);
 o = obj.Hi * tau * (w_in - R*w_ot - g_tilde);
 end
 closure = @closure_fun;
 end
+function closure = boundary_condition_char(obj,boundary,w_data)
-function closure = boundary_condition_outflow(obj, boundary, p_data)
 [e_s,e_S,s] = obj.get_boundary_ops(boundary);
-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)
-% Extract solution at the boundary
 q_s = e_S'*q;
 rho = q_s(1);
 u = q_s(2)/rho;
 e = q_s(3);
 c = obj.c(q_s);
-% Calculate transformation matrix
+Lambda = [u, u+c, u-c];
+p_in = find(s*Lambda < 0);
+p_ot = find(s*Lambda >= 0);
+p = [p_in p_ot];
+pt(p) = 1:length(p);
 T = obj.T(q_s);
-Tin = T(:,p_in);
-Tot = T(:,p_ot);
+tau1 = -2*diag(Lambda(p_in));
+tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char.
-% Convert bc from ordinary form to characteristic form.
-%   Lq = g  => w_in = Rw_ot + g_tilde
-L = (gamma -1)*[0 -1/2*u 1];
-g = [p_data(t)];
-R =-inv(L*Tin)*(L*Tot);
-g_tilde = inv(L*Tin)*g;
-% Setup the penalty parameter
-tau1 = [-2*(c - abs(u))];
-tau2 = [0; 0]; % Penalty only on ingoing char.
 tauHat = [tau1; tau2];
 tau = -s*e_S*sparse(T*tauHat(pt,:));
-% Calculate the numerical caracteristics and apply the penalty
 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);
+w_s_data = w_data(t);
-o = obj.Hi * tau * (w_in - R*w_ot - g_tilde);
+w_in_data = w_s_data(p_in);
+o = obj.Hi * tau * (w_in - w_in_data);
 end
 closure = @closure_fun;
+end
+function closure = boundary_condition_inflow(obj, boundary, p_data, v_data)
+[~,~,s] = obj.get_boundary_ops(boundary);
+switch s
+case -1
+p_in = [1 2];
+case 1
+p_in = [1 3];
+end
+a = obj.gamma - 1;
+L = @(~,u,~)[
+0        1 0;
+0 -1/2*u*a a;
+];
+g = @(t)[
+v_data(t);
+p_data(t);
+];
+closure = boundary_condition_L(obj, boundary, L, g, p_in);
+end
+function closure = boundary_condition_outflow(obj, boundary, p_data)
+[~,~,s] = obj.get_boundary_ops(boundary);
+switch s
+case -1
+p_in = 2;
+case 1
+p_in = 3;
+end
+a = obj.gamma -1;
+L = @(~,u,~)a*[0 -1/2*u 1];
+g = @(t)[p_data(t)];
+closure = boundary_condition_L(obj, boundary, L, g, p_in);
 end
 % Set wall boundary condition v = 0.
 function closure = boundary_condition_wall(obj,boundary)

Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 60:e707cd9e6d0a