Enforcing forward-differencing for compressible Euler Equations

[taylorhughmorgan]

I was able to implement solving the 1D compressible Euler equation in spherically symmetric coordinates using the following code. The code runs and solves the problem, but because it is using central differencing, it quickly becomes unstable. Is there a way to force forward-differencing scheme on all spatial discretization?

System of equations:
$\frac{\partial \rho}{\partial t} + u \cdot \nabla \rho + \rho \nabla \cdot u = 0$
$\frac{\partial u}{\partial t} + u \nabla \cdot u + \frac{\nabla p}{\rho} = 0$
$\frac{\partial e}{\partial t} + u \cdot \nabla e + \frac{p}{\rho} \nabla \cdot u = 0$

code:

from pde import FieldCollection, PDEBase, SphericalSymGrid, ScalarField, VectorField, MemoryStorage, plot_kymograph

class CompressibleEuler(PDEBase):
    """ Solve the compressible Euler equation in spherical symmetry """
    
    def __init__(self, grid, 
                 gamma=1.4, 
                 cv__JpkgK=718, 
                 R__JpkgK =287,
                 method='forward',
                 bc={'derivative' : 0}):
        super().__init__()
        self.bc = bc
        self.gamma = gamma
        self.cv    = cv__JpkgK
        self.R     = R__JpkgK
        self.grid  = grid
        #self.diffOperator = operators.make_derivative(grid, method=method)
        
    def get_initial_state(self, v0, rho0, e0):
        """ prepare a useful, initial state """
        ## converting pressure to phi
        
        v   = VectorField(self.grid, v0, label='velocity')
        rho = ScalarField(self.grid, rho0, label='density')
        e   = ScalarField(self.grid, e0, label='internal energy')
        
        return FieldCollection([v, rho, e])
        
    def evolution_rate(self, state, t=0):
        #assert state.grid.dim == 1  # ensure the state is one-dimensional
        v, rho, e = state
        
        grad_v   = v.gradient(bc=self.bc)
        grad_rho = rho.gradient(bc=self.bc)
        grad_e   = e.gradient(bc=self.bc)
        
        ## relate internal energy back to pressure assuming constant volume specific heat
        p = rho * self.R * e / self.cv
        grad_p = p.gradient(bc=self.bc)
        pOverRho = p / rho
        
        v_t   = -v.dot(grad_v) - grad_p / rho
        rho_t = -v.dot(grad_rho) - rho * v.divergence(bc=self.bc) #rho.dot(grad_v)
        e_t   = -v.dot(grad_e) - pOverRho * v.divergence(bc=self.bc)
        
        return FieldCollection(([v_t, rho_t, e_t]))

if __name__ == '__main__':
    rMin = 0.01
    rMax = 5
    nGridPts = 48
    
    grid  = SphericalSymGrid(radius=[rMin, rMax], shape=nGridPts)  # generate grid
    rGrid = np.linspace(rMin, rMax, num=nGridPts)
    
    ## initial conditions - scaled to size
    u0     = np.zeros_like(rGrid)
    dens0  = np.ones_like(rGrid)
    energ0 = np.ones_like(rGrid)
    energ0[:8] = 3.0
    
    eq = CompressibleEuler(grid)
    storage = MemoryStorage()
    
    state = eq.get_initial_state(u0, dens0, energ0)
    
    sol = eq.solve(state, t_range=1, dt=1e-5, tracker=storage.tracker(0.05))
    sol.plot(kind="image")
    sol.plot()
    plot_kymograph(storage, field_index=2)

[david-zwicker]

Thank you for your message. Forward differences were indeed not very well exposed in the package. I just merged a new pull request in the main branch, which improves the situation. If you download the latest version of the package from GitHub, you should be able to replace your calls to .gradient(bc=self.bc) with .gradient(bc=self.bc, method="forward") to obtain forward derivatives. Analogously, you can use method="backward" to obtain backwards derivatives.

[taylorhughmorgan]

Awesome. Thank you. I did run it, and I can confirm that both method='forward' and method='backward' worked. I have attached files to compare that the 1st-order differencing scheme is indeed more stable.

As a follow up, are there any plans to incorporate this with the divergence operator?
Backwards differencing:

Central differencing:

[david-zwicker]

I just added forward and backward difference schemes for Cartesian and spherical grids.

[taylorhughmorgan]

Awesome. Thank you. I was able to run both the divergence and gradient operators in using the forward, backward, and central differencing methods.