1D hydrostatic equilibrium with Fipy

[NorbertMagyar]

I have originally asked this question on stackoverflow.com, but jeguyer suggested that I move it here. Here is my original description:

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/dx = -\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.

EDIT: The equations that I'm trying to solve are the following (conservation of mass, momentum, and total energy):

$\frac{{\partial}\rho}{\partial t} = - \nabla \cdot (\rho v_x)$
$\frac{{\partial}\rho v_x}{\partial t} = - \nabla \cdot (\rho v_x^2) - \nabla p - \rho g$
$\frac{{\partial}E}{\partial t} = - \nabla \cdot (v_x E + v_x p) - v_x \rho g$

where the pressure $p = (\gamma-1)(E - \rho v_x^2/2)$ is a derived quantity.

Note that in the code written above I have instead implemented the equation for (specific) internal energy:

$\frac{{\partial}E_{th}}{\partial t} = - \nabla \cdot (v_x E_{th}) - (\gamma-1) E_{th} \nabla \cdot v_x$

related to pressure as $p = (\gamma-1)E_{th}$, but meanwhile I believe that the equation for total energy should be implemented.

The boundary conditions would be the following: free flow in/out of the system, thus Neumann-type zero gradient for momentum. Density can also be the same, to allow changes. However, the pressure should be stratified according to the hydrostatic formula in the boundary to prevent draining of the system, at both ends of the 1D domain, therefore its gradient has to be given by $dp/dx = -\rho g$. However, if we do not implement the internal energy as one of the equations but the total energy equation instead, I am not entirely sure how to implement this hydrostatic boundary condition for the total energy.

[guyer]

Can you clarify, is $dp/dz = -\rho g$ both the expected solution and a boundary condition?

[NorbertMagyar]

Sorry, it should be $dp/dx = -\rho g$, as we have $x$-direction. Well, when the system settles, and the velocity is zero, it should be satisfied everywhere. At the boundary it has to be satisfied because otherwise there are forces acting on the boundary due to gravity, thus mass can leak out.

[guyer]

OK. I just wanted to be sure that it should be applied as a boundary condition, thank you.

[NorbertMagyar]

I have also omitted to add the sink term due to gravity to the total energy equation. Corrected now.

[guyer]

Since your equations are hyperbolic, FiPy is going to have trouble.

Flow is not my thing and I've previously had trouble trying to figure out what "outflow" means in the context of our Stokes problem.

I am generally distrustful of 1D formulations; when written as partial derivatives, it's impossible to tell the difference between divergences and gradients, which need to be discretized differently. When written with vector operators, as you have, you end up with divergences of scalars, which makes my head hurt. TLDR, I'd recommend writing your equations in vector form and only reduce to the 1D form at the very end.

[TomVeeDee]

Hi Guyer,

I'm Tom, collaborator of Norbert on this topic.

Could you explain in more detail what you mean with the following statements in your post?

Since your equations are hyperbolic, FiPy is going to have trouble.

Do you mean that fipy is not designed for hyperbolic systems? On the first page of fipy, it is stated that "fipy is an object oriented, partial differential equation (PDE) solver", and it seems there is no subset of the PDEs is preferred. We have also applied fipy to other hyperbolic PDEs that are part of this project, and there fipy seems to work beautifully.
So, could you explain why fipy is going to have trouble with these hyperbolic equations specifically? I would prefer not to abandon fipy, since it seems to be working mostly for other parts of the code (the system of differential equations and the code we have written is much larger than the above minimum working example).

I'd recommend writing your equations in vector form and only reduce to the 1D form at the very end.

Could you also further explain what you mean with this? Do you mean we should write the equations in vector form here? The equations can be written in vector form, and they are equal to the equations on wikipedia that you have mentioned in your answer on stackexchange. For completeness, here they are with the notation we use in the code

$$\frac{\partial \rho}{\partial t} = - \nabla\cdot(\rho \vec{v}),$$
$$\frac{\partial \rho \vec{v}}{\partial t} = - \nabla\cdot(\rho \vec{v}\vec{v}) - \nabla ((\gamma-1) E_{th}) - \rho\vec{g},$$
$$\frac{\partial E_{th}}{\partial t} = -\nabla\cdot(E_{th}\vec{v})-(\gamma-1) E_{th} \nabla\cdot\vec{v},$$
where $\vec{v}\vec{v}$ is the dyadic product of the velocity vector $\vec{v}$. In the code, we use the variable mom to represent the momentum $\rho\vec{v}$. The variable dens is the gas density $\rho$, and the variable Eth is the thermal energy $E_{th}$.

Or, do you mean with your comment that we should write the python code using a vector notation? We can have a try, but could you give a first example how to do it for (e.g.) the continuity equation (eq_continuity)?

Thanks already for your input.

Best regards,

Tom Van Doorsselaere

[guyer]

I'm Tom, collaborator of Norbert on this topic.

Good to meet you, Tom

Could you explain in more detail what you mean with the following statements in your post?

Since your equations are hyperbolic, FiPy is going to have trouble.

Do you mean that fipy is not designed for hyperbolic systems? On the first page of fipy, it is stated that "fipy is an object oriented, partial differential equation (PDE) solver", and it seems there is no subset of the PDEs is preferred. We have also applied fipy to other hyperbolic PDEs that are part of this project, and there fipy seems to work beautifully. So, could you explain why fipy is going to have trouble with these hyperbolic equations specifically? I would prefer not to abandon fipy, since it seems to be working mostly for other parts of the code (the system of differential equations and the code we have written is much larger than the above minimum working example).

FiPy is capable of solving parabolic, elliptical, and hyperbolic PDEs, as well as combinations thereof. I've used it for simple Beer's Law calculations. However, the discretizations we use can mean that it's not terribly accurate (see our advection examples; the examples themselves are a train-wreck and need to be rewritten, but they illustrate the issue). Further, I don't think we can handle shocks at all. We often point people toward Clawpack if they have exclusively hyperbolic equations in demanding parameter spaces. @wd15 is more versed in this stuff and can hopefully pipe up if I'm steering you wrong.

That said, I'm glad that FiPy is working for you more broadly and am not trying to chase you away. I think getting this particular example is working is just a matter of getting the correct formulation for the boundary condition; I think it's probably some sort of Robin condition, but that's why I'm being pedantic about the vector formulation.

I'd recommend writing your equations in vector form and only reduce to the 1D form at the very end.

Could you also further explain what you mean with this? Do you mean we should write the equations in vector form here?
:
Or, do you mean with your comment that we should write the python code using a vector notation? We can have a try, but could you give a first example how to do it for (e.g.) the continuity equation (eq_continuity)?

Sorry for not being clear. I meant the first. If you were doing a multi-dimensional solution, it might be advantageous to write your momentum as a FiPy vector, but it's not required. Since it's essentially a scalar in 1D, I wouldn't do it.

The equations can be written in vector form, and they are equal to the equations on wikipedia that you have mentioned in your answer on stackexchange. For completeness, here they are with the notation we use in the code:

OK, but @NorbertMagyar later said that a total energy statement was better and that it needed a sink term due to gravity. I tried to make the transformation from the internal energy equation on wikipedia, but got myself confused (although I was fighting off a headache at the time). Just so we're all on the same page, can you confirm the energy equation you actually want to use?

[TomVeeDee]

I'm not particularly attached to either energy equation. @NorbertMagyar says the total energy equation is more reliable in case we want to achieve conservation of energy in the system, which is somewhat important. The vector formulation of that equation is

$$ \frac{\partial E}{\partial t} = -\nabla\cdot(\vec{v} (E+p)) - \rho \vec{v}\cdot\vec{g}.$$

I had originally chosen the equation for thermal energy, because of the ease of implementation (our full code also has other energy sources), but we could switch easily.

[TomVeeDee]

Sorry, suddenly I can't get the Latex code to render properly.

[guyer]

The math needed whitespace around it. I'll think about this and get back to you.

[TomVeeDee]

Today, I've gone deeper into this. I've found that the boundary conditions are not the cause of the non-equilibrium, but rather the calculation of the gradient in the momentum equation. In the boundary points, both the gaussGrad or the leastSquaresGrad would have a value which does not match the required value of the boundary condition in the respective position. This then leads to spurious (i.e. numerical) forces, which destabilises the system.

I'm pasting the full code below, but the key point is that I've changed the momentum equation from

eq_momentum   = (fipy.TransientTerm(var=mom) 
            == - fipy.ConvectionTerm(coeff=v*[[1.]],var=mom)
              - (gamma-1)*Eth.leastSquaresGrad[0]
              - fipy.ImplicitSourceTerm(coeff=g,var=dens)
                )

to

tmpterm = ((gamma-1)*Eth.leastSquaresGrad[0]).value
tmpterm[0]= -g*basevalue*np.exp(-mesh.cellCenters[0].value[0]/scaleheight)
tmpterm[-1]= -g*basevalue*np.exp(-mesh.cellCenters[0].value[-1]/scaleheight)
gradp = fipy.CellVariable(name="gradp", mesh=mesh, value=tmpterm, hasOld=True)
eq_momentum   = (fipy.TransientTerm(var=mom) 
            == - fipy.ConvectionTerm(coeff=v*[[1.]],var=mom)
              - gradp 
              - fipy.ImplicitSourceTerm(coeff=g,var=dens)
                )

The latter is a dirty hack, where I replace the first and last point of the gradient with the gradient that is expected in this condition. With that hack, the atmosphere remains stable for around 1000 iterations. I suspect after that time other numerical errors creep in, but that is not the issue that is being discussed here.

Could you suggest a more proper way to include the gradient in the momentum equation? It seems the gradient is numerically computed from the boundary value, and not the boundary gradient. It is this latter that is constrained in this example.

Best regards,

Tom

"""
Smallest workable example of hydrostatic equilibrium. Norbert 22.01.2024
"""

import fipy
import numpy as np
import pylab as plt


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

g = 2.74e4 * unit_length/unit_velocity**2 # test with constant g 

T_ini = 1.e6/unit_temperature # in K 

    
dens=fipy.CellVariable(name="dens",mesh=mesh,value=0., hasOld=True)
scaleheight = 5.e9/unit_length #in cm
scaleheight = T_ini*unit_temperature/(g*unit_velocity**2/unit_length)* kB/mp/a*b /unit_length
basevalue = 1.e-15/unit_density # in g cm^{-3}
dens.setValue(basevalue*np.exp(-mesh.cellCenters/scaleheight))
dens.faceValue.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[0]/scaleheight), where=mesh.facesLeft)
dens.faceValue.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[-1]/scaleheight), where=mesh.facesRight)
        
vr=fipy.CellVariable(name="vr",mesh=mesh,value=0.)
basevalue_v = 0.e0/unit_velocity # cm/s
vr.setValue(basevalue_v*np.exp(vr.mesh.cellCenters/scaleheight))
Ek=dens*pow(vr,2)/2 
    
mom = fipy.CellVariable(name="mom", mesh=mesh, value=0., hasOld=True)
mom.setValue(dens*vr)
mom0=dens.value[0]*vr.value[0]
mom.faceValue.constrain(0.0, where=mesh.facesLeft)
mom.faceValue.constrain(0.0, where=mesh.facesRight)

v=(mom/dens).arithmeticFaceValue


divv=v.divergence

Eth = fipy.CellVariable(name="Eth", mesh=mesh, value=dens.value*T_ini/(gamma-1), hasOld=True)

E0=Eth.value[0]
Eth.faceGrad.constrain(-g*basevalue*np.exp(-mesh.faceCenters[0].value[0]/scaleheight)/(gamma-1), where=mesh.facesLeft)
Eth.faceGrad.constrain(-g*basevalue*np.exp(-mesh.faceCenters[0].value[-1]/scaleheight)/(gamma-1), where=mesh.facesRight)
      
eq_continuity = (fipy.TransientTerm(var=dens) 
            == - fipy.ConvectionTerm(coeff=v*[[1.]], var=dens)
                ) 
eq_eth        = (fipy.TransientTerm(var=Eth) 
            == - fipy.ConvectionTerm(coeff=v*[[1.]],var=Eth) 
               - fipy.ImplicitSourceTerm(coeff=divv,var=Eth)*(gamma-1)
                )
tmpterm = ((gamma-1)*Eth.leastSquaresGrad[0]).value
tmpterm[0]= -g*basevalue*np.exp(-mesh.cellCenters[0].value[0]/scaleheight)
tmpterm[-1]= -g*basevalue*np.exp(-mesh.cellCenters[0].value[-1]/scaleheight)
gradp = fipy.CellVariable(name="gradp", mesh=mesh, value=tmpterm, hasOld=True)
eq_momentum   = (fipy.TransientTerm(var=mom) 
            == - fipy.ConvectionTerm(coeff=v*[[1.]],var=mom)
              # - (gamma-1)*Eth.leastSquaresGrad[0]
              - gradp 
              - fipy.ImplicitSourceTerm(coeff=g,var=dens)
                )
    
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.5, datamax=1.,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=5000

# 60% of maximal sound speed in the domain
cflcond = 0.6 * dx/np.max(np.sqrt(gamma**2*Eth.value/dens.value))
timeStepDuration = cflcond
t=np.arange(steps)*timeStepDuration
    
for step in np.arange(steps):
        Eth.updateOld()
        dens.updateOld()
        mom.updateOld()
        gradp.updateOld()
        eqn.solve(dt=timeStepDuration)

        if step%100 ==0:
            print('done step',step, ', time', step*timeStepDuration*unit_time)
        viewer.plot()

[guyer]

The following changes retain your solution and are clearer (to me, anyway):

diff --git a/euler.py b/euler.py
index 523af86..8fff740 100644
--- a/euler.py
+++ b/euler.py
@@ -46,8 +46,8 @@ scaleheight = 5.e9/unit_length #in cm
 scaleheight = T_ini*unit_temperature/(g*unit_velocity**2/unit_length)* kB/mp/a*b /unit_length
 basevalue = 1.e-15/unit_density # in g cm^{-3}
 dens.setValue(basevalue*np.exp(-mesh.cellCenters/scaleheight))
-dens.faceValue.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[0]/scaleheight), where=mesh.facesLeft)
-dens.faceValue.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[-1]/scaleheight), where=mesh.facesRight)
+dens.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[0]/scaleheight), where=mesh.facesLeft)
+dens.constrain(basevalue*np.exp(-mesh.faceCenters[0].value[-1]/scaleheight), where=mesh.facesRight)
         
 vr=fipy.CellVariable(name="vr",mesh=mesh,value=0.)
 basevalue_v = 0.e0/unit_velocity # cm/s
@@ -57,8 +57,8 @@ Ek=dens*pow(vr,2)/2
 mom = fipy.CellVariable(name="mom", mesh=mesh, value=0., hasOld=True)
 mom.setValue(dens*vr)
 mom0=dens.value[0]*vr.value[0]
-mom.faceValue.constrain(0.0, where=mesh.facesLeft)
-mom.faceValue.constrain(0.0, where=mesh.facesRight)
+mom.constrain(0.0, where=mesh.facesLeft)
+mom.constrain(0.0, where=mesh.facesRight)
 
 v=(mom/dens).arithmeticFaceValue
 
@@ -68,8 +68,7 @@ divv=v.divergence
 Eth = fipy.CellVariable(name="Eth", mesh=mesh, value=dens.value*T_ini/(gamma-1), hasOld=True)
 
 E0=Eth.value[0]
-Eth.faceGrad.constrain(-g*basevalue*np.exp(-mesh.faceCenters[0].value[0]/scaleheight)/(gamma-1), where=mesh.facesLeft)
-Eth.faceGrad.constrain(-g*basevalue*np.exp(-mesh.faceCenters[0].value[-1]/scaleheight)/(gamma-1), where=mesh.facesRight)
+Eth.faceGrad.constrain(-g*dens.faceValue / (gamma-1), where=mesh.exteriorFaces)
       
 eq_continuity = (fipy.TransientTerm(var=dens) 
             == - fipy.ConvectionTerm(coeff=v*[[1.]], var=dens)
@@ -78,14 +77,9 @@ eq_eth        = (fipy.TransientTerm(var=Eth)
             == - fipy.ConvectionTerm(coeff=v*[[1.]],var=Eth) 
                - fipy.ImplicitSourceTerm(coeff=divv,var=Eth)*(gamma-1)
                 )
-tmpterm = ((gamma-1)*Eth.leastSquaresGrad[0]).value
-tmpterm[0]= -g*basevalue*np.exp(-mesh.cellCenters[0].value[0]/scaleheight)
-tmpterm[-1]= -g*basevalue*np.exp(-mesh.cellCenters[0].value[-1]/scaleheight)
-gradp = fipy.CellVariable(name="gradp", mesh=mesh, value=tmpterm, hasOld=True)
 eq_momentum   = (fipy.TransientTerm(var=mom) 
             == - fipy.ConvectionTerm(coeff=v*[[1.]],var=mom)
-              # - (gamma-1)*Eth.leastSquaresGrad[0]
-              - gradp 
+              - (gamma-1)*Eth.leastSquaresGrad[0]
               - fipy.ImplicitSourceTerm(coeff=g,var=dens)
                 )
     
@@ -107,7 +101,6 @@ for step in np.arange(steps):
         Eth.updateOld()
         dens.updateOld()
         mom.updateOld()
-        gradp.updateOld()
         eqn.solve(dt=timeStepDuration)
 
         if step%100 ==0:

The first two sets of changes are just a reflection that when you write var.constrain(value, where=face_locations), FiPy automatically interprets that as var.faceValue.constrain(value, where=face_locations) and it saves a bit of typing.

The next and substantive change is to express the constrain as closely to the math as possible.

$$\begin{align} \left.\frac{d p}{d x}\right|_{\text{inlet}|\text{outlet}} &= -\rho g \\\ \left(\gamma - 1\right)\left.\frac{E_{th}}{d x}\right|_{\text{inlet}|\text{outlet}} &= -\rho g \end{align}$$

or

$$\texttt{Eth.faceGrad.constrain}( -\rho g / \left(\gamma - 1\right), \texttt{where=}\text{inlet}|\text{outlet})$$

You don't need to calculate the value at the specific points because you've already constrained rho to have the appropriate value at those points.

At that point, you don't need gradp for anything.

[TomVeeDee]

Your suggestion does not seem to fix the problem I'm experiencing. I've changed the Eth variable to now have the boundary condition

Eth.faceGrad.constrain(-g*dens.faceValue / (gamma-1), where=mesh.exteriorFaces)

as you suggested, but the gradient of Eth still contains non-physical values in the first and last point. I'm attaching a figure that shows the difference between my hacked variable gradp (which gives the desired physical result) and the Eth.leastSquaresGrad, after just one iteration of the code. It is the unphysical jump of the latter in the first and last point that lead to the unphysical result.

From my side, it seems there is a mismatch between the implementation of the grad methods and the imposed faceGrad boundary condition.

[guyer]

I suspect that this notebook I put together several years ago is relevant: interpolationsDiscretizationsAndGradientsOhMy.ipynb

I had not noticed that you're using leastSquaresGrad. We never use it and I don't know why we have it. As my notebook shows, it doesn't see constraints because of the way it works. I would use Eth.grad, but I see that the boundary velocity is then ill-behaved.

[TomVeeDee]

I think you are right in saying that the current representation of the system is probably the culprit. What is often done to achieve that is to reformulate the

$$ \nabla p $$

as

$$ \nabla \cdot (p \mathbf{I}), $$

where $\mathbf{I}$ is the unit tensor of the correct dimension. Such a formulation would bring it closer to the fipy.ConvectionTerm
However, coding up this term as

- (gamma-1)*fipy.ConvectionTerm(coeff=[[1.]],var=Eth)

does not seem to be the correct understanding on my side.

Could you refer me to an example implementation of a ConvectionTerm written in terms of the unit tensor?

[guyer]

(gamma-1)*fipy.CentralDifferenceConvectionTerm(coeff=[[1.]],var=Eth) seems to produce stable solutions. Eth.grad (and Eth. leastSquaresGrad) still have boundary issues, as illustrated in the notebook I linked above, but they're not used in the solution of the equations.

[TomVeeDee]

Indeed, with this term, it seems stabler, but the solution is not stationary in my case. When running the code with the CentralDifferenceConvectionTerm, the interior of the system is behaving well, i.e. the velocity remains at 0. This shows that this is probably close to the implementation I need for the $\nabla\cdot p\vec{I}$ term.
However, at the left and right edges of the system, I see a systematic decrease of the velocity. Physically speaking, that should not be the case. The velocity (red line in the graph) should remain 0 everywhere at all time.
As a reference, I'm pasting again the python implementation here.

eq_momentum   = (fipy.TransientTerm(var=mom) 
            == - fipy.ConvectionTerm(coeff=v*[[1.]],var=mom)
              - (gamma-1)*fipy.CentralDifferenceConvectionTerm(coeff=[[1]],var=Eth)
              - fipy.ImplicitSourceTerm(coeff=g,var=dens)
                )

Is there a way to fix the 'boundary conditions' for the CentralDifferenceConvectionTerm? I guess it is similarly influenced as the Eth.grad[0] term?

I've tried another implementation (inspired by your python notebook on showing the grad boundary problems).

R = mesh.facesRight
L = mesh.facesLeft
LR = mesh.facesRight | mesh.facesLeft
n = mesh.faceNormals
cellsNearRight = (R * n).divergence
cellsNearLeft = (L * n).divergence
cellsNearEdge = (LR * n).divergence
cellsNearEdge/=np.max(cellsNearEdge.value)
cellsNearRight/=np.max(cellsNearRight.value)
cellsNearLeft/=np.max(cellsNearLeft.value)
eq_momentum   = (fipy.TransientTerm(var=mom) 
            == - fipy.ConvectionTerm(coeff=v*[[1.]],var=mom)
              - (gamma-1)*(Eth.grad[0]-Eth.grad[0]*cellsNearEdge) + g*dens*cellsNearLeft+ g*dens*cellsNearRight
              - fipy.ImplicitSourceTerm(coeff=g,var=dens)
                )

However, this implementation becomes numerically unstable around 300 time steps. So, that hacked solution is also not the correct way to go.

[guyer]

I haven't forgotten about you!

Using examples.flow.stokesCavity and examples.reactiveWetting.liquidVapor1D as guides, I reformulated to use a pressure (energy) corrector and Rhie-Chow interpolation, both of which we know to be necessary for other flow problems. I've posted my notebook attempting to solve this. The systematic decrease in velocity is no longer present, but velocity ringing sets in immediately, and eventually propagates to the other fields.

Note that I shifted from solving for momentum to solving for velocity, as I believe this is necessary for applying the pressure corrector to the continuity equation. I'd have thought the Rhie-Chow interpolation would have dealt with the ringing, but I obviously don't know what I'm doing. Smaller time steps, more sweeps, relaxation, and different types of ConvectionTerm don't help. Maybe @wd15 has some insights.