1D hydrostatic equilibrium with Fipy

Question

I am trying to implement the compressible Euler equations including a full energy equation in Fipy. In the code below, I aim to obtain a hydrostatic solution (dp/dz = -rho g) in the presence of gravity. However, the code does not seem to do the job properly. I am not sure what could be wrong with these boundary conditions:

import fipy
import numpy as np
import pylab as plt

from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

#system of units (relative to CGS): 
# user specified units
unit_length = 1.e8 # cm
unit_numberdensity = 1.e9 # cm^-3
unit_temperature = 1.e6 # K

# constants
He_abundance = 0.1 
mp = 1.67e-24 # g
kB =  1.3806488e-16 # erg K^-1
miu0 = 4*np.pi
gamma=5./3.
R_sun = 6.957e10/unit_length # cm
a = 1.0 + 4.0*He_abundance
b = 2.0 + 3.0*He_abundance
g = 2.74e4 * pow(R_sun,2)/(pow(x+R_sun,2))*unit_length/unit_velocity**2; g.name='g'
T_ini = 1.e6/unit_temperature # in K 
# derived units
unit_density = a*mp*unit_numberdensity 
unit_pressure = b*unit_numberdensity*kB*unit_temperature
unit_velocity = (unit_pressure/unit_density)**0.5
unit_time = unit_length/unit_velocity
scaleheight = (T_ini/g.value) #in cm
dens0 = 1.e-15/unit_density # in g cm^{-3}

L = 1.e10/unit_length # in cm
nx = 1000
dx = L/nx
mesh = fipy.Grid1D(dx=L / nx, nx=nx)

x = mesh.cellCenters[0]
x.name='x'
    
dens=fipy.CellVariable(mesh=mesh,value=dens0*np.exp(-mesh.cellCenters[0]/scaleheight), hasOld=True)
dens.faceValue.constrain(dens0, mesh.facesLeft)
dens.faceValue.constrain(dens.value[-1], mesh.facesRight)

mom = fipy.CellVariable(mesh=mesh, value=0., hasOld=True)
mom.faceGrad.constrain(0, mesh.facesLeft)
mom.faceGrad.constrain(0, mesh.facesRight)

gradpL = g.value[0]*dens.value[0]/(gamma-1)
gradpR = g.value[-1]*dens.value[-1]/(gamma-1)

Eth = fipy.CellVariable(mesh=mesh, value=dens.value*T_ini/(gamma-1), hasOld=True)
Eth.faceGrad.constrain(gradpL, where=mesh.facesLeft)
Eth.faceGrad.constrain(gradpR, where=mesh.facesRight)
      
divv=(mom/dens).arithmeticFaceValue.divergence

eq_continuity = (fipy.TransientTerm(var=dens) 
            == - fipy.ConvectionTerm(coeff=mom*[[1.]]/dens, var=dens)
                ) 
eq_momentum   = (fipy.TransientTerm(var=mom) 
            == - fipy.ConvectionTerm(coeff=mom*[[1.]]/dens,var=mom)
              - (gamma-1)*Eth.grad.dot([[1.]])
              - fipy.ImplicitSourceTerm(coeff=g,var=dens)
                )
eq_eth        = (fipy.TransientTerm(var=Eth) 
            == - fipy.ConvectionTerm(coeff=mom*[[1.]]/dens,var=Eth) 
               - fipy.ImplicitSourceTerm(coeff=divv,var=Eth)*(gamma-1)
                )

eqn = eq_continuity & eq_momentum & eq_eth
     
viewer = fipy.viewers.matplotlibViewer.matplotlib1DViewer.Matplotlib1DViewer(vars=(dens, mom/dens*unit_velocity/1.e6, Eth*(gamma-1)),
                    datamin=-0.2, datamax=0.6,legend='upper left', title = '1D hydro solution')
viewer.axes.legend([r'Density ($10^{-12}\ \mathrm{kg\ m^{-3}}$)',r'Velocity ($\mathrm{10\ km/s}$)',r'Pressure (0.03 Pa)'])
viewer.plot()

steps=10000

# 60% of maximal sound speed in the domain
cflcond = 0.2 * dx/np.max(np.sqrt(gamma**2*Eth.value/dens.value))
timeStepDuration = cflcond
t=np.arange(steps)*timeStepDuration
    
solution=np.zeros((3,steps,len(x.globalValue)))

for step in np.arange(steps):
        Eth.updateOld()
        dens.updateOld()
        mom.updateOld()
        eqn.solve(dt=timeStepDuration)

        if rank==0:
            if step%50 ==0:
                print('done step',step, ', time', step*timeStepDuration*unit_time)
                viewer.plot()
            solution[0,step,:]=dens.globalValue
            solution[1,step,:]=(mom/dens).globalValue
            solution[2,step,:]=Eth.globalValue
        
            if (step%500==0) & (step != 0):
                np.savez_compressed('savefile'+str(step).zfill(int(np.ceil(np.log10(steps+1)))), solution[:,step-500:step,:])

I am expecting the 1D system to reach equilibrium, which does not seem to happen.