[mathgra]

Hello Im trying to write a matlab code for 1d Euler Equation I have tried different codding types. Lastly I have showed 2D code and I have tried to write it in 1d with Local Lax Friedrich method but I thing my result does not converge cause it doesnt show in figure.

Can u help me to find the mistake in my code ?

gamma = 1.4;
CFL = 0.5;
tMax = 0.75;
nx = 11;
dx = 1/(nx-1);

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

Q = zeros(nx,3);
Qnew = zeros(nx,3);
F = zeros(nx,3);
Fh = zeros(nx-1,3);
c = zeros(nx,1);
ah1 = zeros(nx,1);

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

for i=1:nx
if (x(i) < 0.0 )
p(i) = 1.0;
r(i) = 1.0;
u(i) = 0.0;
elseif (x(i) > 0.0 )
p(i) = 0.1;
r(i) = 0.125;
u(i) = 0.0;
end
c(i) = sqrt((gamma*p(i))/r(i));
end

Eig1 = abs(u) + c;
dt = CFL*dx./max(Eig1(: );
nt = floor(tMax/dt);

for k=1:nt
for i=1:nx

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

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=1:nx-1

ah1(i) = max(abs(u(i)) + c(i), abs(u(i+1)) + c(i+1));

Fh(i,: ) = 0.5*(F(i+1,: ) + F(i,: ) - ah1(i).*(Q(i+1,: ) - Q(i,: )));

end

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


for i=2:nx-1
Qnew(i,= Q(i,: ) - (dt/dx)*(Fh(i,: ) - Fh(i-1,: ));
end


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

end

[selig5576]

You will have to be more specific, but this code looks like the I posted some time ago . Your getting NaN solution correct? Try placing the boundary conditions AFTER you have solver for the interior and that should do the trick.

EDIT: Also tMax should be 0.15 or so for Sod's shock tube.

EDIT 2: Also the way you calculate dt is not complete as you still need more conditions. For now just use a fixed time step when debugging.

[mathgra]

Originally Posted by selig5576

You will have to be more specific, but this code looks like the I posted some time ago . Your getting NaN solution correct? Try placing the boundary conditions AFTER you have solver for the interior and that should do the trick.

EDIT: Also tMax should be 0.15 or so for Sod's shock tube.

EDIT 2: Also the way you calculate dt is not complete as you still need more conditions. For now just use a fixed time step when debugging.

Yes it is really similar to your code also. Only boundary conditions different. I have studied your codes first I think you use Python to write code. It take some time for me to understand and then I have found similar code for 2D here and I have mixed two and try to write this.

Yes I get Nan solution I will try your advices thank u so much

[moshe]

This is just personal suggestion, but I think in terms of learning you should attempt things on your own. For me it has generated me to ask many questions, but I think its good to first try it on your own.

[mathgra]

Originally Posted by moshe

This is just personal suggestion, but I think in terms of learning you should attempt things on your own. For me it has generated me to ask many questions, but I think its good to first try it on your own. On a side note, why are you mixing codes? My python code 100% works.

Im still learning Im really beginner so Im studying each part of the code and try to understand how matlab understand the code when we write.
I also try to write different type of codes for this problem. You are right, I agree everything you write. I dont know how to use Python so I use matlab for now.

[rewol]

Well, i solved this problem couple of years ago with FVM, utilizing ROE Flux splitting methodology.

Reviewing the code, I seem to choose my dt as 0.001 with tmax = 1.2.

I remember getting many NaN values in my solution due to my incorrect formulation of eigenvalues and eigenvectors. In MATLAB, usually the reason for NaN is because a dominator in formulation is assigned or calculated as 0 somehow, which was my case anyway.

I usually add isreal function in between the calculations to see where it goes sideways. Because i remember solution vector becomes complex too with NaN values.

Regards.

Alper

[crazyshock]

gamma = 1.4;
CFL = 0.5;
tMax = 0.2;
nx = 101;
dx = 1/(nx-1) ;

r = zeros(nx,1) ;
u = zeros(nx,1) ;
p = zeros(nx,1) ;

Q = zeros(nx,3) ;
F = zeros(nx,3) ;

x=linspace(0,1,nx) ;
half=uint8(nx/2);
p(1:half)=1.0; r(1:half)=1.0; u(: )=0.0;
p(half+1:end)=0.1; r(half+1:end)=0.125;
c=sqrt(gamma*p./r);

Eig1 = abs(u) + c;
dt = CFL*dx/max(Eig1(: ));
nt = floor(tMax/dt);

for k=1:nt

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

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

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


ah1=max(abs(u+c),abs(circshift(u,-1)+circshift(c,-1)));
ah2=repmat(ah1,[1,3]);
Fh = 0.5*( F+circshift(F,-1) - ah2.*(circshift(Q,-1)-Q) );

Qnew = Q - (dt/dx)*(Fh-circshift(Fh,1));


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


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

plot(x,r);
pause(0.1);

end




cheers. Here circshift is used to provide periodic boundary conditions. Also, learn the difference between .* and *

[javeria]

hello can you post your codes

[javeria]

i need ur codes can u post it please

[poly_tec]

na, go so the work yourself.

[javeria]

please I m stuck

[poly_tec]

Where are you stuck at?

[javeria]

I used laxwendroff and macckormack schemes , approximately reimann solver..... but now I need some other scheme any high resolution scheme to complete my research ... so I m stuck now