Mercurial > repos > public > sbplib

comparison +scheme/Euler1d.m @ 50:75ebf5d3cfe5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Performance improvement for Euler1d
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 05 Nov 2015 17:10:57 -0800
parents 8f0c2dc747dd
children 2194cd385419
comparison
equal deleted inserted replaced
49:8f0c2dc747dd 50:75ebf5d3cfe5
12 M % Derivative norm 12 M % Derivative norm
13 gamma 13 gamma
14 14
15 T 15 T
16 p 16 p
17 c
17 18
18 19
19 H % Discrete norm 20 H % Discrete norm
20 Hi 21 Hi
21 e_l, e_r, e_L, e_R; 22 e_l, e_r, e_L, e_R;
77 % q=[rho, rho*u, e]^T 78 % q=[rho, rho*u, e]^T
78 % F=[rho*u, rho*u^2+p, (e+p)*u] ^T 79 % F=[rho*u, rho*u^2+p, (e+p)*u] ^T
79 % p=(gamma-1)*(e-rho*u^2/2); 80 % p=(gamma-1)*(e-rho*u^2/2);
80 81
81 %Solving on form q_t + F_x = 0 82 %Solving on form q_t + F_x = 0
82 function o = F(q) 83 function o = F(Q)
83 o = [q(2); q(2).^2/q(1) + p(q); (q(3)+p(q))*q(2)/q(1)]; 84 % Flux: f = [q2; q2.^2/q1 + p(q); (q3+p(q))*q2/q1];
85 o = [Q(2,:); Q(2,:).^2/Q(1,:) + p(Q); (Q(3,:)+p(Q)).*Q(2,:)./Q(1,:)];
84 end 86 end
85 87
86 % Equation of state 88 % Equation of state
87 function o = p(q) 89 function o = p(Q)
88 o = (gamma-1)*(q(3)-q(2).^2/q(1)/2); 90 % Pressure p = (gamma-1)*(q3-q2.^2/q1/2)
89 end 91 o = (gamma-1)*(Q(3,:)-Q(2,:).^2/Q(1,:)/2);
92 end
93
94 function o = c(Q)
95 % Speed of light c = sqrt(obj.gamma*p/rho);
96 o = sqrt(gamma*p(Q)./Q(1,:));
97 end
98
90 99
91 % R = 100 % R =
92 % [sqrt(2*(gamma-1))*rho , rho , rho ; 101 % [sqrt(2*(gamma-1))*rho , rho , rho ;
93 % sqrt(2*(gamma-1))*rho*u , rho*(u+c) , rho*(u-c) ; 102 % sqrt(2*(gamma-1))*rho*u , rho*(u+c) , rho*(u-c) ;
94 % 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]); 103 % 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]);
108 ]; 117 ];
109 % Devide columns by stuff to get rid of extra comp? 118 % Devide columns by stuff to get rid of extra comp?
110 end 119 end
111 120
112 function o = Fx(q) 121 function o = Fx(q)
113 o = zeros(size(q)); 122 Q = reshape(q,3,m);
114 for i = 1:3:3*m 123 o = reshape(F(Q),3*m,1);
115 o(i:i+2) = F(q(i:i+2));
116 end
117 o = D1*o; 124 o = D1*o;
118 end 125 end
119 126
120 127
121 % A=R*Lambda*inv(R), där Lambda=diag(u, u+c, u-c) (c är ljudhastigheten) 128 % A=R*Lambda*inv(R), där Lambda=diag(u, u+c, u-c) (c är ljudhastigheten)
131 obj.Fx = @Fx; 138 obj.Fx = @Fx;
132 obj.T = @T; 139 obj.T = @T;
133 obj.u = x; 140 obj.u = x;
134 obj.x = kr(x,ones(3,1)); 141 obj.x = kr(x,ones(3,1));
135 obj.p = @p; 142 obj.p = @p;
143 obj.c = @c;
136 obj.gamma = gamma; 144 obj.gamma = gamma;
137 end 145 end
138 146
139 147
140 % Enforces the boundary conditions 148 % Enforces the boundary conditions
252 260
253 p = [p_in, p_zero, p_ot]; % Permutation to sort 261 p = [p_in, p_zero, p_ot]; % Permutation to sort
254 pt(p) = 1:length(p); % Inverse permutation 262 pt(p) = 1:length(p); % Inverse permutation
255 263
256 function o = closure_fun(q) 264 function o = closure_fun(q)
257 p = obj.p(q);
258 265
259 q_s = e_S'*q; 266 q_s = e_S'*q;
260 rho = q_s(1); 267 rho = q_s(1);
261 u = q_s(2)/rho; 268 u = q_s(2)/rho;
262 c = sqrt(obj.gamma*p/rho); 269 c = obj.c(q_s);
263 270
264 T = obj.T(q_s); 271 T = obj.T(q_s);
265 R = -(u-c)/(u+c); 272 R = -(u-c)/(u+c);
266 % l = [u, u+c, u-c]; 273 % l = [u, u+c, u-c];
267 274