[mathgra]

Hello Im working on numerical solution of Sod Shock tube problem. I have written matlab code using Lax Friedrichs. Now I wanna write matlab code with MacCormack two step scheme.

My code is working but atfer some steps I got NaN values so I can not obtain my figures. If you can look my code and give any idea for solving this problem I will be thankful. And density start at 0.5 instead of 1 I dont understand why this happen.

Here the code that I write in MATLAB ;

clear all;
close all;
clc;

nx = 101;

% Initialization
r = zeros(nx,1);
u = zeros(nx,1);
E = zeros(nx,1);
p = zeros(nx,1);

x = zeros(nx,1);
Q = zeros(nx,3);
F = zeros(nx,3);
Qnew = zeros(nx,3);
c = zeros(nx,1);

% Parameters
gamma = 1.4;
CFL = 0.68;
tMax = 0.2;
dx = 1/(nx-1);

% Creating mesh

for i=1:nx
x(i) = (i-1)*dx;
end


% Initial condition
for i=1:nx
if (x(i) > 0.5)
p(i) = 0.1;
r(i) = 0.125;
u(i) = 0.0;

elseif (x(i) <= 0.5 )
p(i) = 1.0;
r(i) = 1.0;
u(i) = 0.0;

end
E(i) = p(i)./(gamma-1) + 0.5*r(i)*(u(i).^2);
c(i) = sqrt(gamma*p(i)/r(i));
end


% Q_new will be the solution vector

Qnew(:,1)=r;
Qnew(:,2)=r.*u;
Qnew(:,3)=E;

t=0.d0;
% Time loop
while (t<tMax)
dt = CFL*dx./(max(max(c+sqrt(u.*u))));
for i=1:nx

% Flux calculation
Q(i,1) = r(i);
Q(i,2) = r(i).*u(i);
Q(i,3) = E(i);


F(i,1) = r(i).*u(i);
F(i,2) = r(i).*u(i).^2 + p(i);
F(i,3) = u(i).*(E(i) + p(i));
end




for i=2:nx-1

% predictor step
Qstar(i,1) = Q(i,1) - (dt/dx)*(F(i+1,1) - F(i,1)) ;
Qstar(i,2) = Q(i,2) - (dt/dx)*(F(i+1,2) - F(i,2)) ;
Qstar(i,3) = Q(i,3) - (dt/dx)*(F(i+1,3) - F(i,3)) ;

% Corrector step

ustar(i) = Qstar(i,2) / Qstar(i,1);
pstar(i) = (gamma-1)*( Qstar(i,3) - 0.5*(Qstar(i,2).^2/Qstar(i,1)) );

Fstar(i,1) = Qstar(i,1).*ustar(i);
Fstar(i,2) = Qstar(i,1).*(ustar(i).^2) + pstar(i);
Fstar(i,3) = ustar(i).*( Qstar(i,3) + pstar(i) );

Qnew(i,1) = 0.5* (Qnew(i,1) + Qstar(i,1) - (dt/dx) * (Fstar(i,1) - Fstar(i-1,1)));

Qnew(i,2) = 0.5* (Qnew(i,2) + Qstar(i,2) - (dt/dx) * (Fstar(i,2) - Fstar(i-1,2)));

Qnew(i,3) = 0.5* (Qnew(i,3) + Qstar(i,3) - (dt/dx) * (Fstar(i,3) - Fstar(i-1,3)));


end


% Boundary condition

Qnew(1,: ) = Qnew(2,: );
Qnew(nx,: ) = Qnew(nx-1,: );

% Primitive variables
r = Qnew(:,1);
u = Qnew(:,2)./r;
E = Qnew(:,3);
p = (gamma-1)*(E-0.5*r.*(u.*u));

t=t+dt; % Advance in time

end

plot(x,r);
figure(1)
xlabel('x');
title('Density');

[naffrancois]

I think I can see one coding mistake. In the corrector step you use Fstar(i-1,: ). As the components of Fstar are defined in a loop from 2 to nx, Fstar(1,: ) is never computed.

Second comment, probably not critical here, in general it is better to update the speed of sound inside the temporal loop for dt computation

[mathgra]

Originally Posted by naffrancois

I think I can see one coding mistake. In the corrector step you use Fstar(i-1,: ). As the components of Fstar are defined in a loop from 2 to nx, Fstar(1,: ) is never computed.

Second comment, probably not critical here, in general it is better to update the speed of sound inside the temporal loop for dt computation

Oh you are right, Thanks a lot I added Fstar at boundary as

Fstar(1,: ) = Fstar(2,: );
Fstar(nx,: ) = Fstar(nx-1,: );

but density still start with 0.8 its really interesting

[mathgra]

Originally Posted by naffrancois

I think I can see one coding mistake. In the corrector step you use Fstar(i-1,: ). As the components of Fstar are defined in a loop from 2 to nx, Fstar(1,: ) is never computed.

Second comment, probably not critical here, in general it is better to update the speed of sound inside the temporal loop for dt computation

I did it

thank you so much for advise the place that I put Fstar boundary was wrong.

Now it works Im so much thankful

[naffrancois]

All right well done!

[volkanyildirim48]

Hi mathgra,


I am trying to use your code to have a basic understanding of the mccormick solution, but i keep having error "Undefined function or variable 'Fstar'.
"


I don't know what I am doing wrong.


Thank you

[black key]

Please share your updated code for Maccormack. I am also working on the same for understanding.