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?


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

    os.system('rm results/*')

Lx,Ly = (30,30)

RHO = 12
DT = 0.02 # Main integration timestep
H = 1

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

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

# Enforce the Courant-Friedrichs-Lewy condition
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]))
    # 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)
    #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)

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')

  • $\begingroup$ 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. $\endgroup$ Jan 16 at 9:01

1 Answer 1


There is a lot of one can mention here, so I am going to kick off things.

  1. CFL is never a sufficient criterion for stability, but a necessary one. That said, just because you adhere the CFL condition, there is no guarantee that you get stable results. This makes absolute sense, since your methods can be arbitrarily bad.

  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}$$

I am not really sure what you are doing with the RK1 quantities, i.e., how you want to obtain the gradients. Maybe it's best to start with Euler for the time discretization and (central) finite differences for the spatial derivatives to get a first, simple version.

  • $\begingroup$ Thanks Dan, I re-wrote it to solve it with the Euler method and indeed it yielded acceptable solutions without tending towards becoming unstable. Thank you for the clarity, could you provide more info (ref?) for the claim "the Euler equations are intrinsically for the compressible case"? Thanks. $\endgroup$
    – Gaston
    Jan 16 at 11:41
  • $\begingroup$ From Wikipedia: However, fluid dynamics literature often refers to the full set – including the energy equation – of the more general compressible equations together as "the Euler equations" $\endgroup$
    – Dan Doe
    Jan 16 at 11:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.