# Direct integration of 2D Euler Equations with Runge Kutta shows oscillating Courant-Friedrichs-Lewy coefficient. Stiff or Bug?

By writing the direct integration of the 2D Euler Equations in a wide and short box where the fluid enters and exits through the horizontal faces using the Runge Kutta O(4) method I have found that the CFL coefficient, here approximately computed as follows, oscillates with an increasing amplitude over time. If integration is not stopped, the velocity field diverges.

Am I facing the situation commonly referred to as a direct solver being unstable for a stiff equation?

$$CFL_{coefficient}=\frac{max(u,v)dt}{h}$$

Where $$u,v$$ are the speed in x and y respectively

$$h$$ the spacial characteristic length

$$dt$$ the integration step.

The density is taken as constant.

A couple of notes. I tried enforcing zero tangential speed in the borders, as a humble approximation to simulating viscosity without explicitly adding the viscous term, and it appeared to promote the instability. It wasn't clear tho.

The main note would be considering that the conservation of mass, as follows, is not being enforced at all.

$$\frac{\partial u }{\partial x} + \frac{\partial v }{\partial y} = 0$$

A minimum working example in Python3 is as follows. Please do not run it in any space where a folder called "results" already exists. Requirements are numpy and matplotlib.

EDIT: Indeed the answer is that there was a bug. I simplified the solver to use the O(1) Euler method as suggested by the accepted answer. I have updated the MWE. The timestep is now computed as to satisfy the CLF coeff at every integration step.

import numpy as np
import matplotlib.pyplot as plt
import os,sys

try:
os.mkdir('results')
except:
pass
try:
os.system('rm results/*')
except:
pass

Lx,Ly = (30,30)

RHO = 12
DT = 0.02 # Main integration timestep
H = 1
RHO_H = RHO / H
BOUNDARYVAL = 0.1

Ux = np.ones((Lx, Ly)) * BOUNDARYVAL / 2
Uy = np.zeros((Lx, Ly))

Ux[0,:] = BOUNDARYVAL
#Ux[1,:] = BOUNDARYVAL/2
#Ux[2,:] = BOUNDARYVAL/4

# Enforce the Courant-Friedrichs-Lewy condition
CFL = BOUNDARYVAL * DT / H
assert(CFL <= 0.5)

LAPS = 300
FREQ = 15
D = len(str(LAPS))
for _ in range(LAPS*FREQ):

if _>0:
DT = CFL * H / max([max(Uy.flatten()), max(Ux.flatten())])

# Down to up velocity vertically
Ux[0,:] = BOUNDARYVAL
#Ux[-1,:] = BOUNDARYVAL

# Vertical velocity in the borders can be whatever because there is no viscosity
#Ux[:,-1] = 0
#Ux[:,0] = 0

# Horizonal velocity in the horizontal borders can be whatever
#Uy[0,:] = 0
#Uy[-1,:] = 0

# Horizontal velocity in the vertical borders must be zero as "fluid cant escape through the sides"
Uy[:,0] = 0
Uy[:,-1] = 0

# Plot
if _%FREQ==0:
print(f"DT is {DT}")
print(f"Lap {_+1} of {FREQ*LAPS}!")
plt.imshow(Ux * Ux + Uy * Uy, cmap = 'coolwarm', vmin=min([BOUNDARYVAL,0]), vmax=max([BOUNDARYVAL,0]))
plt.colorbar()
plt.savefig(f'results/{str(_//FREQ).zfill(D)}.png')
plt.close()

# Compute next state
# Mass conservation eq. is divergence(U)=0
# Momentum conservation eqs. are dU/dt + (U \cdot \grad) U = 0
# Where there is no gravity, viscosity, nor external pressure gradient being enforced

# dU/dt = -RHO * (Ux d/dx + Uy d/dy)U
#
# so we compute auxiliaries
#
dUx_dx = (Ux[2:,1:-1] - Ux[:-2,1:-1]) / 2 / H
dUx_dy = (Ux[1:-1,2:] - Ux[1:-1,:-2]) / 2 / H
dUy_dx = (Uy[2:,1:-1] - Uy[:-2,1:-1]) / 2 / H
dUy_dy = (Uy[1:-1,2:] - Uy[1:-1,:-2]) / 2 / H
#
# and evolve the field
#
#print(dUx_dx.shape, Ux.shape)
#plt.imshow(dUx_dx);plt.show()
#MX = np.max((DT * RHO * (Ux[1:-1,1:-1] * dUx_dx + Uy[1:-1,1:-1] * dUx_dy)).flatten())
#print(f'max incoming update is {MX}')
Ux[1:-1,1:-1] -= DT * RHO * (Ux[1:-1,1:-1] * dUx_dx + Uy[1:-1,1:-1] * dUx_dy)
Uy[1:-1,1:-1] -= DT * RHO * (Ux[1:-1,1:-1] * dUy_dx + Uy[1:-1,1:-1] * dUy_dy)

try:
os.remove('result.avi')
except: pass

print(f'To make a GIF showing the evolution please run ffmpeg  -framerate 5 -i "results/%0{D}d.png" -vcodec mpeg4 result.avi')


• This is not the Runge-Kutta method. You need to repeat the steps you did to compute RK1 from the current state also for RK2,... starting from the intermediate states. That is, if RK1 = F(U), then RK2 = F(U+DT/2*RK1) etc. This can be decomposed into 2 steps, construction of the intermediate state U_tmp=U+DT/2*RK1 and the evaluation RK2=F(U_tmp). You only did the first, but even that not correctly, and somehow assigned its result as the result of the second step. Jan 16 at 9:01

2. Looks to me you are trying not to solve the "typical" Euler Equations (which are intrinsically for the compressible case) $$\partial_t \begin{pmatrix} \rho \\ \rho \boldsymbol{u} \\ E \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \boldsymbol{u} \\ \rho \boldsymbol{u}^T \boldsymbol{u} \\ (E + p) \boldsymbol{u} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}$$ but instead something like incompressible Navier-Stokes with neglected viscosity, pressures and body forces: $$\nabla \cdot \boldsymbol{u} = 0, \quad \partial_t \boldsymbol{u} + (\boldsymbol{u} \cdot \nabla) \boldsymbol{u} = \boldsymbol{0}$$ The formulation above is overdetermined, since we have $$1 +2$$ equations, but only two unknowns: $$u_x, u_y$$. So okay, just as you did let's not enforce conservation of mass but just evolve $$\boldsymbol{u}$$ according to the momentum equation: \begin{align} u_x &= - u_x\partial_x u_x - u_y\partial_y u_x \\ u_y &= -u_x\partial_x u_y - u_y\partial_y u_y \end{align}