Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 61:183d9349b4c1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactored euler1d to have real methods.
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 14 Nov 2015 19:26:30 -0800
parents e707cd9e6d0a
children 00344139cd7d
comparison
equal deleted inserted replaced
60:e707cd9e6d0a 61:183d9349b4c1
9 9
10 D % non-stabalized scheme operator 10 D % non-stabalized scheme operator
11 M % Derivative norm 11 M % Derivative norm
12 gamma 12 gamma
13 13
14 T
15 p
16 c
17 Lambda
18
19
20 H % Discrete norm 14 H % Discrete norm
21 Hi 15 Hi
22 e_l, e_r, e_L, e_R; 16 e_l, e_r, e_L, e_R;
23 17
24 end 18 end
74 % Man har Q_t+F_x=0 i 1D Euler, där 68 % Man har Q_t+F_x=0 i 1D Euler, där
75 % q=[rho, rho*u, e]^T 69 % q=[rho, rho*u, e]^T
76 % F=[rho*u, rho*u^2+p, (e+p)*u] ^T 70 % F=[rho*u, rho*u^2+p, (e+p)*u] ^T
77 % p=(gamma-1)*(e-rho*u^2/2); 71 % p=(gamma-1)*(e-rho*u^2/2);
78 72
73
79 %Solving on form q_t + F_x = 0 74 %Solving on form q_t + F_x = 0
80 function o = F(Q)
81 % Flux: f = [q2; q2.^2/q1 + p(q); (q3+p(q))*q2/q1];
82 o = [Q(2,:); Q(2,:).^2./Q(1,:) + p(Q); (Q(3,:)+p(Q)).*Q(2,:)./Q(1,:)];
83 end
84
85 % Equation of state
86 function o = p(Q)
87 % Pressure p = (gamma-1)*(q3-q2.^2/q1/2)
88 o = (gamma-1)*( Q(3,:)-1/2*Q(2,:).^2./Q(1,:) );
89 end
90 obj.p = @p;
91
92 function o = c(Q)
93 % Speed of light c = sqrt(obj.gamma*p/rho);
94 o = sqrt(gamma*p(Q)./Q(1,:));
95 end
96 obj.c = @c;
97
98 function o = Lambda(q)
99 u = q(2)/q(1);
100 c = obj.c(q);
101 L = [u, u+c, u-c];
102 o = diag(L);
103 end
104 obj.Lambda = @Lambda;
105
106 function o = T(q)
107 % T is the transformation such that A = T*Lambda*inv(T)
108 % where Lambda=diag(u, u+c, u-c)
109 rho = q(1);
110 u = q(2)/q(1);
111 e = q(3);
112
113 c = sqrt(gamma*p(q)/rho);
114
115 sqrt2gamm = sqrt(2*(gamma-1));
116
117 o = [
118 sqrt2gamm*rho , rho , rho ;
119 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ;
120 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 ;
121 ];
122 % Devide columns by stuff to get rid of extra comp?
123 end
124 obj.T = @T;
125 75
126 function o = Fx(q) 76 function o = Fx(q)
127 Q = reshape(q,3,m); 77 Q = reshape(q,3,m);
128 o = reshape(F(Q),3*m,1); 78 o = reshape(obj.F(Q),3*m,1);
129 o = D1*o; 79 o = D1*o;
130 end 80 end
131 81
132 function o = Fx_disp(q) 82 function o = Fx_disp(q)
133 Q = reshape(q,3,m); 83 Q = reshape(q,3,m);
134 f = reshape(F(Q),3*m,1); 84 f = reshape(obj.F(Q),3*m,1);
135 85
136 c = c(Q); 86 c = obj.c(Q);
137 lambda_max = c+abs(Q(2,:)./Q(1,:)); 87 lambda_max = c+abs(Q(2,:)./Q(1,:));
138 % lambda_max = max(lambda_max); 88 % lambda_max = max(lambda_max);
139 89
140 lamb_Q(1,:) = lambda_max.*Q(1,:); 90 lamb_Q(1,:) = lambda_max.*Q(1,:);
141 lamb_Q(2,:) = lambda_max.*Q(2,:); 91 lamb_Q(2,:) = lambda_max.*Q(2,:);
155 obj.u = x; 105 obj.u = x;
156 obj.x = kr(x,ones(3,1)); 106 obj.x = kr(x,ones(3,1));
157 obj.gamma = gamma; 107 obj.gamma = gamma;
158 end 108 end
159 109
110 % Flux function
111 function o = F(obj, Q)
112 % Flux: f = [q2; q2.^2/q1 + p(q); (q3+p(q))*q2/q1];
113 p = obj.p(Q);
114 o = [Q(2,:); Q(2,:).^2./Q(1,:) + p; (Q(3,:)+p).*Q(2,:)./Q(1,:)];
115 end
116
117 % Equation of state
118 function o = p(obj, Q)
119 % Pressure p = (gamma-1)*(q3-q2.^2/q1/2)
120 o = (obj.gamma-1)*( Q(3,:)-1/2*Q(2,:).^2./Q(1,:) );
121 end
122
123 % Speed of sound
124 function o = c(obj, Q)
125 % Speed of light c = sqrt(obj.gamma*p/rho);
126 o = sqrt(obj.gamma*obj.p(Q)./Q(1,:));
127 end
128
129 % Eigen value matrix
130 function o = Lambda(q)
131 u = q(2)/q(1);
132 c = obj.c(q);
133 L = [u, u+c, u-c];
134 o = diag(L);
135 end
136
137 % Diagonalization transformation
138 function o = T(obj, q)
139 % T is the transformation such that A = T*Lambda*inv(T)
140 % where Lambda=diag(u, u+c, u-c)
141 rho = q(1);
142 u = q(2)/q(1);
143 e = q(3);
144 gamma = obj.gamma;
145
146 c = sqrt(gamma*obj.p(q)/rho);
147
148 sqrt2gamm = sqrt(2*(gamma-1));
149
150 o = [
151 sqrt2gamm*rho , rho , rho ;
152 sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ;
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 ;
154 ];
155 % Devide columns by stuff to get rid of extra comp?
156 end
160 157
161 % Enforces the boundary conditions 158 % Enforces the boundary conditions
162 % w+ = R*w- + g(t) 159 % w+ = R*w- + g(t)
163 function closure = boundary_condition(obj,boundary, type, varargin) 160 function closure = boundary_condition(obj,boundary, type, varargin)
164 [e_s,e_S,s] = obj.get_boundary_ops(boundary); 161 [e_s,e_S,s] = obj.get_boundary_ops(boundary);