# Difficult bug in my 2D Compressible Euler solver

For the past few days, I have been writing a numerical solver for the 2D compressible Euler equations for an ideal gas. My numerical method has been the Local Lax Friedrichs or "Rusanov's method." Before writing the 2D code I successfully wrote a 1D code using the Local Lax-Friedrichs method, which works. However, my extension to 2D has given me some problems I am struggling with no success. In terms of the specific error, I keep getting NaNs in my arrays (for $$\rho$$,$$u$$, $$v$$, $$E$$, $$p$$). I have attached my code below.

To better understand whats going on in the code the method I am using is

\begin{align} Q_{i,j}^{n+1} &= Q_{i,j}^{n} - \frac{\Delta t}{\Delta x} \left(F_{i+1/2,j}^{n} - F_{i-1/2,j}^{n} \right) - \frac{\Delta t}{\Delta y} \left(G_{i,j+1/2}^{n} - G_{i,j-1/2}^{n}\right)\\ F_{i+1/2,j}^{n} &= \frac{1}{2} \left[ \left(F_{i+1,j}^{n} + F_{i,j}^{n}\right) - a_{i+1/2,j} \left(Q_{i+1,j}^{n} - Q_{i,j}^{n} \right)\right]\\ a_{i+1/2,j} &= \max{\left(\left|u_{i,j}\right| + c_{i,j}, \left|u_{i+1,j}\right| + c_{i+1,j}\right)}\\ G_{i,j+1/2}^{n} &= \frac{1}{2} \left[ \left(G_{i,j+1}^{n} - G_{i,j}^{n}\right) - b_{i,j+1/2}\left(Q_{i,j+1}^{n} - Q_{i,j}^{n}\right) \right]\\ b_{i,j+1/2} &= \max{\left(\left|v_{i,j}\right| + c_{i,j}, \left|v_{i,j+1}\right| + c_{i,j+1}\right)} \end{align}

In terms of $$c_{i,j}$$, this term refers to the speed of sound, i.e. $$c_{i,j} = \sqrt{\frac{\gamma p}{\rho}}$$ Our pressure can be defined as

$$p = (\gamma - 1)(E - \frac{\rho}{2} \left(u^{2} + v^{2}\right)\\ E = \frac{p}{\gamma-1} + \frac{\rho}{2} \left(u^{2} + v^{2}\right)\\ Q = [\rho\;\rho u\;\rho v\;E]^{T}\\ F = [\rho u\; \rho u^{2} + p\; \rho u v\; u(E+p)]^{T}\\ G = [\rho v\; \rho u v\; \rho v^{2} + p\; v(E+p)]^{T}$$

nx = 101;
dx = 1/(nx-1);

ny = 101;
dy = 1/(ny-1);

r = zeros(nx,ny);
u = zeros(nx,ny);
v = zeros(nx,ny);
E = zeros(nx,ny);
p = zeros(nx,ny);

x = zeros(nx,ny);
y = zeros(nx,ny);

Q = zeros(nx,ny,4);
Qnew = zeros(nx,ny,4);

F = zeros(nx,ny,4);
G = zeros(nx,ny,4);

Fh = zeros(nx-1,ny,4);
Gh = zeros(nx,ny-1,4);

gamma = 1.4;

c = zeros(nx,ny);
ah1 = zeros(nx,ny);
ah2 = zeros(nx,ny);

CFL = 0.5;
tMax = 0.25;

for j=1:ny
for i=1:nx
x(i,j) = (i-1)*dx;
y(i,j) = (j-1)*dy;
end
end

for j=1:ny
for i=1:nx
if (x(i,j) >= 0.5 && y(i,j) >= 0.5)
p(i,j) = 0.4;
r(i,j) = 0.5313;
u(i,j) = 0.0;
v(i,j) = 0.0;
elseif (x(i,j) < 0.5 && y(i,j) >= 0.5)
p(i,j) = 1.0;
r(i,j) = 1.0;
u(i,j) = 0.7276;
v(i,j) = 0.0;
elseif (x(i,j) < 0.5 && y(i,j) < 0.5)
p(i,j) = 1.0;
r(i,j) = 0.8;
u(i,j) = 0.0;
v(i,j) = 0.0;
elseif (x(i,j) >= 0.5 && y(i,j) < 0.5)
p(i,j) = 1.0;
r(i,j) = 1.0;
u(i,j) = 0.0;
v(i,j) = 0.7276;
end

c(i,j) = sqrt((gamma*p(i,j))/r(i,j));
end
end

Eig1 = abs(u) + c;
Eig2 = abs(v) + c;

dt = CFL*dx./max(Eig1(:),Eig2(:));
nt = floor(tMax/dt);

for k=1:nt
for j=1:ny
for i=1:nx
E(i,j) = p(i,j)./(gamma-1) + (0.5*r(i,j))*(u(i,j).^2 + v(i,j).^2);

Q(i,j,1) = r(i,j);
Q(i,j,2) = r(i,j).*u(i,j);
Q(i,j,3) = r(i,j).*v(i,j);
Q(i,j,4) = E(i,j);

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

G(i,j,1) = r(i,j).*v(i,j);
G(i,j,2) = r(i,j).*v(i,j).*u(i,j);
G(i,j,3) = r(i,j).*v(i,j).^2 + p(i,j);
G(i,j,4) = v(i,j).*(E(i,j) + p(i,j));
end
end

for j=1:ny-1
for i=1:nx-1
ah1(i,j) = max(abs(u(i,j)) + c(i,j), abs(u(i+1,j)) + c(i+1,j));
ah2(i,j) = max(abs(v(i,j)) + c(i,j), abs(v(i,j+1)) + c(i,j+1));

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

end
end

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

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

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

figure(1)
contour(x,y,r);
xlabel('x');
ylabel('y');
title('Density');


In terms of the initial condition, it divides the discrete domain into 4 regions, i.e. 4 quadrants.

My code does "run," but it only plots the initial condition which means there's something going wrong in my time update (however for my 1D problem it works with no problem.) My density plot should look like this:

What we should expect in terms of ($$\rho$$, $$u$$, $$v$$, $$p$$) is something similar to the picture I uploaded. The problem I am solving is a test problem which tests the accuracy of numerical solvers (Wendroff & Liska 2003)

• Your code works fine for me (MATLAB 2015a). Add clear all at the beginning to be sure to clear all variables when you re-run your code. I indeed get 4 quadrants but can you explain what we should find in term of (r,u,v,p) ? I would have written dt = CFL*min(dx,dy)./max(Eig1(:),Eig2(:)); if the number of nodes in x and y direction is not the same. – Coriolis Aug 14 '16 at 10:17

First, your time step dt is a matrix. That is not possible. Just below your calculate the number of required time step to reach tMax. You do not know that. You do not know how many time step you will require to complete your simulation, it depends on the velocity field and the speed of sound. The time step is linked to the mesh size and the maximum velocity through the CFL condition such that, for the Euler problem : $$dt = CFL \frac{h}{\vert \vert \mathbf{u} \vert \vert + c}$$ where h is the size of the cell. Notice here that we deal with the norm of the velocity, not the absolute value. In 1D, the absolute value works, you were just lucky. To find the largest admissible time step, you need to compute the maximum velocity (that is to say $max(\vert \vert \mathbf{u} \vert \vert + c)$) through your computational domain because each cell will have a different velocity. If your mesh is uniformly spaced, it is the cell with the maximum velocity that will dictate the maximum time step allowed. This process absolutely need to be inside your loop in time at the beginning of each iteration, your dt will vary through time and it is a scalar, not a matrix that will change at each iteration.

Second, your boundary conditions are not properly set. Actually, you iterate through all your cells except ones on the boundaries which have a value of zero at the first iteration because of your initialization of $Q_{new}$. Add keyboard just before the end of the time loop and check your workspace, you will see that u and v are NaN because you have divided by 0.

BC are difficult to handle, a common solution is to add ghost cells with a virtual state so you don't need to separate the computation of the boundary cells, or directly impose the flux on the boundary.

EDIT :

The following piece of code works for me, however I can't ensure you it is bug-free. I usually work with an isothermal Euler model so I don't know if my energy equation is implemented in the right way. Double-check the pressure law and notation between $E$ and $e$ to be sure.

clear all;
close all;
clc;

nx = 101;
ny = 101;

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

x = zeros(nx,ny);
y = zeros(nx,ny);

Q = zeros(nx,ny,4);
F = zeros(nx,ny,4);
G = zeros(nx,ny,4);

Fh = zeros(nx-1,ny,4);
Gh = zeros(nx,ny-1,4);

Qnew = zeros(nx,ny,4);

c = zeros(nx,ny);
ah1 = zeros(nx,ny);
ah2 = zeros(nx,ny);

% Parameters
gamma = 1.4;
CFL = 0.5;
tMax = 0.25;
dx = 1/(nx-1);
dy = 1/(ny-1);

% Creating mesh
for j=1:ny
for i=1:nx
x(i,j) = (i-1)*dx;
y(i,j) = (j-1)*dy;
end
end

% Initial condition
for j=1:ny
for i=1:nx
if (x(i,j) >= 0.5 && y(i,j) >= 0.5)
p(i,j) = 0.4;
r(i,j) = 0.5313;
u(i,j) = 0.0;
v(i,j) = 0.0;
elseif (x(i,j) < 0.5 && y(i,j) >= 0.5)
p(i,j) = 1.0;
r(i,j) = 1.0;
u(i,j) = 0.7276;
v(i,j) = 0.0;
elseif (x(i,j) < 0.5 && y(i,j) < 0.5)
p(i,j) = 1.0;
r(i,j) = 0.8;
u(i,j) = 0.0;
v(i,j) = 0.0;
elseif (x(i,j) >= 0.5 && y(i,j) < 0.5)
p(i,j) = 1.0;
r(i,j) = 1.0;
u(i,j) = 0.0;
v(i,j) = 0.7276;
end
E(i,j) = p(i,j)./(gamma-1) + 0.5*r(i,j)*(u(i,j).^2+v(i,j).^2);
c(i,j) = sqrt(gamma*p(i,j)/r(i,j));
end
end

% Q_new will be the solution vector
Qnew(:,:,1)=r;
Qnew(:,:,2)=r.*u;
Qnew(:,:,3)=r.*v;
Qnew(:,:,4)=E;

t=0.d0;
% Time loop
while (t<tMax)
dt = CFL*min(dx,dy)./(max(max(c+sqrt(u.*u+v.*v)))); % Compute time step
for j=1:ny
for i=1:nx

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

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

G(i,j,1) = r(i,j).*v(i,j);
G(i,j,2) = r(i,j).*v(i,j).*u(i,j);
G(i,j,3) = r(i,j).*v(i,j).^2 + p(i,j);
G(i,j,4) = v(i,j).*(E(i,j) + p(i,j));
end
end

% Rusanov method
for j=1:ny-1
for i=1:nx-1
ah1(i,j) = max(abs(u(i,j)) + c(i,j), abs(u(i+1,j)) + c(i+1,j));
ah2(i,j) = max(abs(v(i,j)) + c(i,j), abs(v(i,j+1)) + c(i,j+1));

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

end
end

% Advance the solution in time
for j=2:ny-1
for i=2:nx-1
Qnew(i,j,:) = Qnew(i,j,:) - (dt/dx)*(Fh(i,j,:) - Fh(i-1,j,:)) - (dt/dy)*(Gh(i,j,:) - Gh(i,j-1,:));
end
end

% Boundary condition
Qnew(:,1,:)  = Qnew(:,2,:);
Qnew(1,:,:)  = Qnew(2,:,:);
Qnew(:,ny,:) = Qnew(:,ny-1,:);
Qnew(nx,:,:) = Qnew(nx-1,:,:);

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

end

% Plot
figure(1)
contour(x(2:nx-1,2:nx-1),y(2:nx-1,2:nx-1),r(2:nx-1,2:nx-1));
xlabel('x');
ylabel('y');
title('Density');


The main changes are the while time loop and the computation of dt. Also I have slightly changed your boundary conditions (moving after the update in time). Actually, your mesh now incorporates ghost cells for cells indexes i=1,i=nx,j=1 and j=ny. These cells are not part of the computational mesh (the update of the solution occurs for 2:nx-1 and 2:ny-1) but you need them to compute flux. I just copy the updated value from the neighbouring cells as you do before. That is not perfect but I can't investigate much further, you didn't mention the boundary condition for your problem (Dirichlet, Neumann ?). That's why in the contour command, I don't plot these cells.

• Hi thank you for the informative answer. However, there are some things I am not sure about. 1. Is there any particular norm I should be using (L1, L2, etc?), 2. In terms of calculating my dt considering my problem is 2D, I am assuming my dt would look like dt = CFL*h/norm(u+c,v+c) since we have u and v velocity. – Simon Aug 16 '16 at 16:07
• Sorry my notation is confusing. $\mathbf{u}$ in bold means a vector so $\mathbf{u}=(u,v)$. The norm is the euclidean norm: $\vert \vert \mathbf{u} \vert \vert =\sqrt[]{u^2+v^2}$ – Coriolis Aug 16 '16 at 16:30
• Hi Coriolis, I've read what you said and have implemented the suggestion of how to calculate dt, but am getting no luck. Aside from the boundary conditions, is there something I am not handling properly? In terms of the calculation of dt, I've tried very small values (that would be within the criteria of convergence) and am getting NaNs, while when I use your suggested dt I don't get any advancement in time. Sorry if I am missing someting – Simon Aug 16 '16 at 17:04
• I fixed your boundary conditions (no NaN anymore) but I still got negative pressure. Could you edit your post to add your closure pressure law ? Mention also the components of the vector Q (primitive or conservative variables ?). Have you validated your 1D implementation with a Sod shock tube ? – Coriolis Aug 16 '16 at 18:46
• Hi Coriolis, Thanks for pointing that out, however from a few places I've looked at, the equations I have are not incorrect (people.nas.nasa.gov/~pulliam/Classes/New_notes/euler_notes.pdf - Page 2). With that said I will try and change it. Would you mind posting your code so I can see if its converging to the correct solution (taking into account the initial conditions and boundary conditions?) – Simon Aug 17 '16 at 13:40