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

\[\tilde{p}_i = p_i - p_{\mathrm{hse},i}\]

This is designed to work with HSE discretized as:

\[p_{i+1} = p_i + \frac{1}{4} \Delta x (\rho_i + \rho_{i+1}) (g_i + g_{i+1})\]

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:

\[\Delta p_i = \Delta \tilde{p}_i + \rho_i g_i \Delta x\]

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 to 0.

  • 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 to 1, 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 requires hse_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).