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

Question

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

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.