0
$\begingroup$

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

enter image description here

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))
ax = fig.add_subplot(gs[0])
cbar_ax = fig.add_subplot(gs[1])

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

$\endgroup$
3
  • $\begingroup$ @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. $\endgroup$
    – F-B
    Sep 25, 2023 at 23:06
  • $\begingroup$ Why would you implement a method from a 50-year old paper? Aren't there better ones available these days? $\endgroup$ Sep 25, 2023 at 23:24
  • $\begingroup$ @WolfgangBangerth Hi, I'm replicating the results of a paper (agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/92JA02905 ) that uses this algorithm. $\endgroup$
    – F-B
    Sep 25, 2023 at 23:56

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.