Where has my implementation of the SHAS2D (Book and Boris, 1975) model gone wrong?

I am new to computational fluid dynamics, and am having a tough time implementing the Book and Boris flux limiter SHAS2D model.

I am attempting to re-create the SHAS2D fluid code to solve the Euler equations. The algorithm for the code is found in the appendix of this paper (Equations A6 to A17):

D.L. Book; J.P. Boris; K. Hain (1975). Flux-corrected transport II: Generalizations of the method. , 18(3), 248–283. doi:10.1016/0021-9991(75)90002-9

https://www.sciencedirect.com/science/article/abs/pii/0021999175900029

To test if I've implemented the algorithm correctly, instead of working with the entire set of Euler equations, I chose to model only the mass continuity equation ($$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0$$) for an arbitrary density with zero velocity everywhere. If $$\vec{v} = 0$$, then $$\frac{\partial \rho}{\partial t} = 0$$ and I'd expect whatever density profile was initially set to stay exactly the same over time.

The set up is summarized as follows:

A 2D box spanning X = [-1,1] and Y = [-1,1], with dx = 2./100., dy = 2./100., and dt = dx/3.

• Density ($$\rho$$) is set to $$\rho$$ = $$e^{-10(x^2 + y^2)}$$
• X - component of momentum density is set to $$\rho \cdot V_x$$ = 0
• Y - component of momentum density is set to $$\rho \cdot V_y$$ = 0

With this set up I expect the density profile to stay the same, however, when I run my code I get the weird result shown below (the plot shows $$\rho (x, y = 0, t = 0)$$ in red and $$\rho (x, y = 0, t = 1000*dt)$$ in blue):

I've checked my code many times, but I can't find where the error is. Where has my implementation gone wrong? Here is my code below:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.gridspec import GridSpec

def SHAST_2D(F,v,S,dl,dt,coord,nx,ny):

Q_p = np.zeros((ny,nx))
Q_n = np.zeros((ny,nx))
if coord == 'x':

F_d = np.zeros((ny,nx))
# step 1 - computing untransported fluxes and the diffused (but untransported) solution
k = np.diff(F,axis = 1) # (eq. A6)
F_d[:,1:-1] = F[:,1:-1] + np.diff(k,axis = 1) # (eq. A7)

# step 2 - computing the transport convection factors
eps_p = 1/2. + v*(dt/dl) # (eq. A8)
eps_n = 1/2. - v*(dt/dl) # (eq. A8)
Q_p[:,:-1] = eps_n[:,:-1]/ (eps_p[:,1:] + eps_n[:,:-1]) # (eq. A9)
Q_n[:,1:] = 1 - Q_p[:,:-1] # (eq. A10)

# step 3 - computing the tranported and diffused solution and changes due to transport
F_td = np.zeros((ny,nx))
F_td[:,1:-1] = (4.*Q_p[:,1:-1]**2.)*k[:,1:] - (4.*Q_n[:,1:-1]**2.)*k[:,:-1] + Q_p[:,1:-1]*(F[:,1:-1] - S[:,1:]) + Q_n[:,1:-1]*(F[:,1:-1] - S[:,:-1]) # (eq. A11)
d_F_t = np.zeros((ny,nx))
d_F_t[:,1:-1] = F_td[:,1:-1] - F_d[:,1:-1] # (eq. A12)

# step 4 - Computing the difference in the transported + diffused solution and initial fluxes (eq. A6) are modified
d_F_td = np.diff(F_td,axis = 1)  # (eq. A13)
k_t = k + 1/8.*np.diff(d_F_t,axis = 1) # (eq. A14)

# step 5 - computing anti-diffusion fluxes
s_kt = np.sign(k_t)
mins = np.min(np.stack((s_kt[:,1:-1] * d_F_td[:,:-2], np.abs(k_t[:,1:-1]), s_kt[:,1:-1]*d_F_td[:,2:])),axis = 0)
k_c = np.zeros(np.shape(k))
k_c[:,1:-1] = s_kt[:,1:-1]*np.max(np.stack((np.zeros(np.shape(mins)), mins)),axis = 0) # (eq. A15)

# step 6 - computing the corrected flux
F_new = F
F_new[:,1:-1] = F_td[:,1:-1] - np.diff(k_c,axis = 1) # (eq. A16)

else:

F_d = np.zeros((ny,nx))
# step 1
k = np.diff(F,axis = 0)
F_d[1:-1,:] = F[1:-1,:] + np.diff(k,axis = 0)

# step 2
eps_p = 1/2. + v*(dt/dl)
eps_n = 1/2. - v*(dt/dl)
Q_p[:-1,:] = eps_n[:-1,:]/ (eps_p[1:,:] + eps_n[:-1,:])
Q_n[1:,:] = 1. - Q_p[:-1,:]

# step 3
F_td = np.zeros((ny,nx))
F_td[1:-1,:] = (4.*Q_p[1:-1,:]**2.)*k[1:,:] - (4.*Q_n[1:-1,:]**2.)*k[:-1,:] + Q_p[1:-1,:]*(F[1:-1,:] - S[1:,:]) + Q_n[1:-1,:]*(F[1:-1,:] - S[:-1,:])
d_F_t = np.zeros((ny,nx))
d_F_t[1:-1,:] = F_td[1:-1,:] - F_d[1:-1,:]

# step 4
d_F_td = np.diff(F_td,axis = 0)
k_t = k + 1/8.*np.diff(d_F_t,axis = 0)

# step 5
s_kt = np.sign(k_t)
mins = np.min(np.stack((s_kt[1:-1,:] * d_F_td[:-2,:], np.abs(k_t[1:-1,:]), s_kt[1:-1,:]*d_F_td[2:,:])),axis = 0)
k_c = np.zeros(np.shape(k))
k_c[1:-1,:] = s_kt[1:-1,:]*np.max(np.stack((np.zeros(np.shape(mins)), mins)),axis = 0)

# step 6
F_new = F
F_new[1:-1,:] = F_td[1:-1,:] - np.diff(k_c,axis = 0)

if np.any(np.isnan(F_new)):
raise ValueError('F is nan')
return F_new

nx = 100
ny = 100

dx = 2./nx
dy = 2./ny
dt = dx/3.

x = np.arange(-1,1,step = dx)
y = np.arange(-1,1,step = dy)
X,Y = np.meshgrid(x,y)

''' --- density --- '''
rho = np.exp(-10.*(X**2. + Y**2.))

''' --- x-momentum density --- '''
rho_vx = np.zeros((ny,nx))

''' --- y-momentum density --- '''
rho_vy = np.zeros((ny,nx))

S_x = np.zeros((ny,nx-1))

frames = 1000
rho_all = np.zeros((ny,nx,frames))
for f in range(frames):

print(f)
vx = rho_vx/rho
vy = rho_vy/rho

rho_new = SHAST_2D(rho,vx,S_x,dx,dt,'x',nx,ny)
rho = rho_new
rho_all[:,:,f] = rho_new

''' ------ for animation ------ '''
gs = GridSpec(1, 2, width_ratios = [0.9, 0.05])
fig = plt.figure(figsize = (7,7))

def myPlot(ax,X,Y,Z):
c = ax.plot(X[50,:],Z)
#c = ax.contourf(X,Y,Z,levels = np.linspace(0,3,100))
#fig.colorbar(c,cax = cbar_ax, shrink = 0.5, aspect = 5)

def update_rho(frame,ax,data):
ax.cla()
#Z = data[:,:,frame]
Z = data[50,:,frame]
myPlot(ax,X,Y,Z)

#ani = animation.FuncAnimation(fig=fig, func=update_rho, frames=frames, fargs = (ax,rho_all),interval=30)
#plt.show()
#ani.save(filename="fluid_example.mp4", writer="ffmpeg")
''' ---------------------------- '''

fig2,ax2 = plt.subplots()
plot1 = ax2.plot(x,rho_all[50,:,-1])
plot2 = ax2.plot(x,rho_all[50,:,0],alpha = 0.5,color = 'red')
ax2.set_ylabel(r'$$\rho (x,y = 0)$$')
ax2.set_xlabel('x')
ax2.legend([r'$$\rho (x,y = 0,t = 1000*dt)$$', r'$$\rho (x,y = 0,t = 0)$$'])


• @njuffa hi, the code uses 64-bit floating point numbers which in python I believe are equivalent to double precision. I increased the time iteration from 1000 to 4000 and found no difference in the final result. I did find one interesting feature however. If I use a square wave (for example rho = 1 for |X| >0.5 and rho = 2 for |X| <= 0.5, then there are no errors. The square wave profile is maintained.
– F-B
Sep 25, 2023 at 23:06
• Why would you implement a method from a 50-year old paper? Aren't there better ones available these days? Sep 25, 2023 at 23:24
• @WolfgangBangerth Hi, I'm replicating the results of a paper (agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/92JA02905 ) that uses this algorithm.
– F-B
Sep 25, 2023 at 23:56