Hydrostatic Equilibrium
There are a few places where we take care to maintain HSE. These differ by integration method:
CTU: this does characteristic tracing of the interface states and is used by Strang and simplified-SDC. CTU requires that we predict the interface states at the midpoint in time – this can make enforcing HSE much more difficult.
true-SDC: this uses a method of lines time reconstruction of the interface states – they are only spatially reconstructed, not extrapolated in time.
Reconstruction
Piecewise linear
Piecewise linear reconstruction can be used with both CTU and true-SDC.
For the 4th-order MC limiter, if we set castro.use_pslope = 1
,
then we subtract off the hydrostatic pressure before limiting.
The idea here is that the hydrostatic pressure balances gravity
so it is not available to generate waves or cause over/undershoots.
We define the excess pressure as \(\tilde{p}\):
This is designed to work with HSE discretized as:
The 4th order MC limiter uses information in zones
\(\{i-2,i-1,i,i+1,i+2\}\). We pick zone i
as the reference
and integrate from i
to the 4 other zone centers, and construct
\(\{\tilde{p}_{i-2}, \tilde{p}_{i-1}, \tilde{p}_{i}, \tilde{p}_{i+1}, \tilde{p}_{i+2}\}\). We then limit on this, giving us the change in the excess
pressure over the zone, \(\Delta \tilde{p}_i\). Finally, we
construct the total slope (including the hydrostatic part) as:
Note
This can be disabled at low densities by setting castro.pslope_cutoff_density
.
The remainder of the PLM algorithm is unchanged. Since this was just the reconstruction part, this step is identical between CTU and true-SDC.
PPM
The PPM version of use_pslope
is essentially the same as the PLM
version described above. We first compute the excess pressure,
\(\tilde{p}\), and then fit a cubic to the 4 cells surrounding an
interface to get the initial interface values. These become the left
and right values of the parabola in each zone. The PPM limiting is
then done on the parabola, again working with \(\tilde{p}\).
Finally, the parabola values are updated to include the hydrostatic
pressure.
We can do better with PPM, and only use the perturbational pressure,
\(\tilde{p}\), in the characteristic tracing and then add back the
hydrostatic pressure to the interface afterwards. This is done via
castro.ppm_well_balanced=1
.
Fully fourth-order method
The 4th order accurate true SDC solver does not appear to need any well-balancing. It maintains HSE to very high precision, as shown in [76].
Boundary conditions
reflecting
For piecewise linear reconstruction, in slope.H
, when using the
4th order MC slopes (castro.plm_limiter = 2
), we impose special
ghost cell values on the normal velocity at reflecting boundaries to
ensure that the velocity goes to 0 at the boundary. This follows
[10] , page 162 (but note that they have a sign
error). Only the normal velocity is treated specially.
HSE
For hydrostatic boundary conditions, we follow the method from [74]. Essentially, this starts with the last zone inside of the boundary and integrates HSE into the ghost cells, keeping the density either constant or extrapolating it. Together the discretized HSE equation and the EOS yield the density and pressure in the ghost cells.
Warning
The HSE boundary condition only works with constant gravity at the moment.
To enable this, we set the appropriate boundary’s xl_ext_bc_type
, xr_ext_bc_type
,
yl_ext_bc_type
, yr_ext_bc_type
, zl_ext_bc_type
,
zr_ext_bc_type
to 1
.
We then control the behavior via the following options.
For the velocity, we have:
hse_zero_vels
: all 3 components of the velocity in the ghost cell are set to0
.hse_reflect_vels
: the normal velocity is reflected and the transverse velocity components are given a zero gradient.
If neither of these are set, then all components of the velocity are simply given a zero gradient.
The temperature in the ghost cells is controlled by:
hse_interp_temp
: if this is set to1
, then we fill the temperatures in the ghost cells via linear extrapolation, using the 2 interior zones just inside the domain. Otherwise, we take the temperature in the ghost cells to be constant.hse_fixed_temp
: if this is positive, then we set the temperature in the ghost cells to the value specified. This requireshse_interp_temp = 0
.
Interface states at reflecting boundary
For all methods, we enforce the reflecting condition on the interface states directly by reflecting the state just inside the domain to overwrite the state on the reflecting boundary just outside of the domain. This is done for all variables (flipping the sign on the normal velocity state). This is especially important for reconstruction that used a one-sided stencil (like the 4th order method).
Test problems
Castro/Exec/gravity_tests/hse_convergence_general
can be used to
test the different HSE approaches. This sets up a 1-d X-ray burst
atmosphere (based on the flame_wave
setup). Richardson
extrapolation can be used to measure the convergence rate (or just
look at how the peak velocity changes).