Compile-Time Flags and Options¶
Most Phlegethon configuration is done at compile time through OPTS += ... entries in each test Makefile.
Makefile options are passed as:
OPTS += USE_MHD
OPTS += nx1_make=64
OPTS += cfl_make=0.8_rp
These are forwarded as preprocessor definitions when compiling source/source.F90.
List of options in Phlegethon¶
Some of the options refer to sections, equations, or tables of the instrument paper.
1. Paths¶
Option |
Meaning |
|---|---|
|
Path to source directory. |
|
Path to |
2. Grid variables¶
Option |
Meaning |
|---|---|
|
Number of spatial dimensions (allowed values are |
|
Global number of grid cells in the |
|
Global number of grid cells in the |
|
Global number of grid cells in the |
|
Number of local domains in the |
|
Number of local domains in the |
|
Number of local domains in the |
|
Number of ghost cells (see Sect. 2.1). |
|
Number of active/passive scalars. |
|
Number of species participating in a nuclear reaction network ( |
|
Number of reactions in the nuclear network. |
2. Grid geometry¶
To set a grid geometry, you need to define one of the following options:
Option |
Meaning |
|---|---|
|
Cartesian geometry (uniform). |
|
Cartesian geometry (nonuniform). |
|
2D polar geometry (requires |
|
2D cylindrical geometry (requires |
|
2D spherical geometry (requires |
|
3D spherical geometry (requires |
|
2D/3D cubed-sphere geometry of Calhoun+08. |
For polar, cylindrical, or spherical grids, the user can enable arbitrarily nonuniform radial grid spacing with NONUNIFORM_RADIAL_NODES.
In this case, the subroutine create_geometry must be defined by the user in app.F90 (see, e.g., tests/logarithmic-grid/app.F90)
The same applies when non-uniform Cartesian grids are enabled with GEOMETRY_CARTESIAN_NONUNIFORM.
For GEOMETRY_2D_CYLINDRICAL, NONUNIFORM_RADIAL_NODES will also enable non-uniform vertical spacing.
For GEOMETRY_CUBED_SPHERE, the following parameters must be defined at compile time:
Option |
Meaning |
|---|---|
|
Maximum radius of the grid. |
|
Regulates the sphericity of the grid. |
|
Regulates the radial stretching of the grid. |
3. Physics modules¶
Option |
Meaning |
|---|---|
|
Solves the governing equations in the co-rotating frame defined by the angular frequency vector |
|
Applies a constant-acceleration body force to the system. In |
|
Enables ideal magnetohydrodynamics (solved with the CT-contact algorithm of Gardiner+05). In |
|
Minimum local Mach number to enable the (upwind) contact version of constrained transport. |
|
Adds magnetic resistivity to the induction equation (note: not to the energy equation and only for Cartesian grids). For this option, the cell-centered value of the magnetic diffusion coefficient |
|
Advects |
|
If declared, advects |
|
Enables space dependent gravity (Newtonian gravity only). The cell-centered components of gravity must be filled in |
|
Sets |
|
Evolves the total energy per unit volume, i.e., \(e \rightarrow e + \phi\). In |
|
Couples a time-dependent nuclear network to the system of conservation laws, in Godunov-split fashion by default (see Sect. 2.13). The network solver is time-implicit. The network has to be used in combination with the option |
|
Enables thermal diffusion solved explicitly in time-unsplit fashion (unstable if the time step is longer than the parabolic CFL criterion, see Sect. 2.10). The cell-centered opacity must be filled in |
|
Enables radiative diffusion solved with the RKL2 super time stepper of Meyer+14 (see, Sect. 2.10). The cell-centered opacity must be filled in |
|
Enables time independent heating. The heating rate per unit volume ( |
|
Switches on a fixed heating source after |
|
Time after which the heating source is activated. |
|
Enables nonnuclear neutrino cooling, computed according to Itoh+1996 (adapted from Frank Timmes’ cococubed). |
4. Boundary conditions (see Table 1)¶
Option |
Meaning |
|---|---|
|
Periodic boundary conditions on the x1 axis. |
|
Periodic boundary conditions on the x2 axis. |
|
Periodic boundary conditions on the x3 axis. |
|
Reflecting boundary conditions at the lower x1 boundary. |
|
Reflecting boundary conditions at the upper x1 boundary. |
|
Reflecting boundary conditions at the lower x2 boundary. |
|
Reflecting boundary conditions at the upper x2 boundary. |
|
Reflecting boundary conditions at the lower x3 boundary. |
|
Reflecting boundary conditions at the upper x3 boundary. |
|
Outflow boundary conditions at the lower x1 boundary. |
|
Outflow boundary conditions at the upper x1 boundary. |
|
Outflow boundary conditions at the lower x2 boundary. |
|
Outflow boundary conditions at the upper x2 boundary. |
|
Outflow boundary conditions at the lower x3 boundary. |
|
Outflow boundary conditions at the upper x3 boundary. |
|
Diode boundary conditions at the lower x1 boundary. |
|
Diode boundary conditions at the upper x1 boundary. |
|
Diode boundary conditions at the lower x2 boundary. |
|
Diode boundary conditions at the upper x2 boundary. |
|
Diode boundary conditions at the lower x3 boundary. |
|
Diode boundary conditions at the upper x3 boundary. |
|
Fixes the temperature (in this example) at the lower x1 boundary only for the thermal diffusion step. If |
|
Imposes reflecting boundary conditions at solid interfaces embedded in the numerical domain. This option only works for Cartesian grids and requires |
With reflecting boundary conditions, the option, e.g., X1L_BFIELD_PMC imposes zero gradient for the normal component and zero for the transverse component of the field.
6. Equation of state (see Sect. 2.2)¶
OPTION A:
By default, Phlegethon uses the EoS of a classical ideal gas with fixed mean molecular weight and adiabatic index, which need to be specified in app.F90.
The option USE_PRAD, adds thermal radiation to the mixture. In app.F90, lgrid%temp(i,j,k) must be filled at every cell center to provide an initial estimate of the temperature for the rhoP->T and rhoe->T transformations appearing in the source code which are based on a Newton–Raphson root finding algorithm.
If ADVECT_SPECIES is enabled, the mean molecular weight is computed from the mass fractions of the species assuming full ionization.
OPTION B:
Alternatively, one can use the Helmholtz EoS (see Timmes+2000):
Option |
Meaning |
|---|---|
|
Uses biquintic interpolation of the free energy to compute the electron-positron EoS contribution (tabulated). In this EoS, electrons/positrons are modeled as a non interacting mixture in an arbitrarily degenerate and relativistic Fermi gas. Additional, ions are assumed to behave like as a classical ideal gas and radiation is in thermal equilibrium. As in the mode |
|
This variable contains the path to the Helmholtz table and needs to be defined every time |
|
Number of nodes in the temperature axis in the Helmholtz table. |
|
Number of nodes in the density axis in the Helmholtz table. |
|
Minimum of \(\mathrm{log}_{10}(T)\) in the Helmholtz table. |
|
Maximum of \(\mathrm{log}_{10}(T)\) in the Helmholtz table. |
|
Minimum of \(\mathrm{log}_{10}(\rho Y_e)\) in the Helmholtz table. |
|
Maximum of \(\mathrm{log}_{10}(\rho Y_e)\) in the Helmholtz table. |
|
Switches on coulomb corrections in the Helmholtz EoS according to Yakovlev+89’s prescription. |
OPTION C:
Alternatively, one can use the PIG Eos of Vetter et al. 2026, in prep.
Option |
Meaning |
|---|---|
|
Computes a partially-ionized-gas EoS using biquintic interpolation of the free energy (PIG, tabulated) + thermal radiation. All ionization states of a bunch of elements and the electron number density are computed using Maxwell–Boltzmann statistics, and radiation is in thermal equilibrium with the gas. As in the mode |
|
This variable contains the path to the PIG EoS table and needs to be defined every time |
|
Number of nodes in the temperature axis in the PIG table. |
|
Number of nodes in the density axis in the PIG table. |
|
Minimum of \(\mathrm{log}_{10}(T)\) in the PIG table. |
|
Maximum of \(\mathrm{log}_{10}(T)\) in the PIG table. |
|
Minimum of \(\mathrm{log}_{10}(\rho)\) in the PIG table. |
|
Maximum of \(\mathrm{log}_{10}(\rho)\) in the PIG table. |
With USE_PRAD, HELMHOLTZ_EOS, or PIG_EOS, the option USE_FASTEOS reconstructs a pair of states for sound speed and internal energy at every grid cell interface to avoid calling the Helmholtz EoS within the Riemann solver subroutine (see Sect. 2.6). If this option is used in combination with USE_WB, the equilibrium state for the thermodynamic variable gammae must be provided at cell centers and cell faces (i.e., lgrid%eq_gammae_cc(i,j,k), lgrid%eq_gammae_x1(i,j,k), etc., see ./tests/hotbubble-2.0-helmholtz).
7. Time Integration (see Sect. 2.9)¶
Option |
Meaning |
|---|---|
|
CFL-based adaptive timestep (Eq. 64). If not used, the code will estimate the time step from |
Choose one of the following time-marching schemes for solving the hyperbolic equations:
Option |
Meaning |
|---|---|
|
Second-order RK time stepper of Shu+88. |
|
Third-order RK time stepper of Shu+88. |
Option |
Meaning |
|---|---|
|
CFL factor |
The option SKIP_HYDRO will skip the hyperbolic update during the simulation.
8. Godunov algorithm¶
Choose one of the following spatial reconstruction methods (see Sect. 2.5)
Option |
Meaning |
|---|---|
|
Second-order linear reconstruction with van Leer slope limiting. |
|
Unlimited parabolic reconstruction for all variables except active scalars, for which limited parabolic reconstruction is used, based on Leidi+24. |
|
Limited fifth-order reconstruction, based on Leidi+24. |
The option USE_SHOCK_FLATTENING will enable the shock-flattening procedure described in Sect. 2.7 (i.e., switch to a two-wave HLL Riemann solver and van Leer reconstruction in the presence of a shock). For this option to be used, the shock-flattening parameter needs to be defined by the user, e.g., eps_sf_make=0.2). This is such that the shock flattener will be activated if the relative jump in pressure at a certain cell location is larger than eps_sf_make.
Riemann solvers (see Sect. 2.6). For hydrodynamic simulations, choose one of the following options:
Option |
Meaning |
|---|---|
|
Three-wave HLLC solver of Toro+94. |
|
Three-wave low-Mach HLLC solver of Minoshima+21. |
For MHD simulations, choose one of the following
Option |
Meaning |
|---|---|
|
Five-wave HLLD solver of Miyoshi+05. |
|
Five-wave low-Mach HLLD solver of Minoshima+21. |
For LHLL-type solvers, a cut-off Mach number must be specified to avoid introducing too-little dissipation in quasi-static flows, e.g., Mach_cutoff_make=1e-4_rp.
For LHLL-type solvers, if both low-Mach and supersonic flows need to be captured on the same grid, the option FLUX_IS_ALL_MACH will switch to a HLL-type solver when Mach>0.6.
9. Gravity options¶
Option |
Meaning |
|---|---|
|
uses the Deviation method of Berberich+19. In |
|
Enables time-dependent gravity by solving Poisson’s equation (see Sect. 2.12), in Godunov-splitting fashion unless |
|
Solves Poisson’s equation for the gravitational potential in every substage of the Runge–Kutta time stepper. |
|
Adds the quadrupole expansion term of the mass distribution to compute boundary conditions for the gravitational potential. |
|
Adds the octupole expansion term of the mass distribution to compute boundary conditions for the gravitational potential (must be used in combination with |
|
Minimum tolerance for the iterative gravity solver, measured in L2 mean error norm (see Eq. 104). |
|
Employs a monopole gravity solver. This option only works for internal boundaries (both 2D and 3D). |
10. Thermal diffusion options¶
Option |
Meaning |
|---|---|
|
Solves the temperature equation rather than the internal energy equation during the radiative diffusion step (see Sect. 2.10). To increase performance, the heat capacity \(c_v\) is held fixed during the time step, so that the EoS does not have to be evaluated in every substep of the super time stepper. |
|
Adaptively changes the number of substeps in the super time stepper. |
|
Balances the thermal diffusion operator using the equilibrium states from the well-balancing method (needs |
|
Enables the computation of radiative+conductive opacities as a function of density, composition and temperature according to Timmes+2000 (implementation basde on coccubed) |
11. Nuclear network options¶
Option |
Meaning |
|---|---|
|
Corrects the reverse rates from reaclib by including the temperature-dependent part of the partition functions of the products and reactants following Rauscher+2000. |
|
Uses the Oda-Langarke-Martinez-Pinedo tabulated weak rates for the nuclear reaction network. |
|
Uses Godunov splitting (1st-order) to couple hydrodynamics to the nuclear network. |
|
Uses Strang splitting (2nd-order) to couple hydrodynamics to the nuclear network. |
|
Uses the backward Euler method (1st-order) to solve the nuclear network. |
|
Uses the esdirk23 method (3 stages, 2nd-order) to solve the nuclear network. |
|
Tolerance of the nonlinear solver used in the nuclear network calculation. |
|
Tolerance of the iterative linear solver (bicgstab) used in the nuclear network calculation. |
|
Minimum temperature to activate the network. |
|
Uses electron screening corrections for reaction rates (2-body reactions only for now), based on Wallace+1982. |
|
Activates boosting of nuclear reactions rates (note: both the rate of change of the abundances and the nuclear energy generation will be boosted). |
|
Boost factor for the nuclear reactions rates. |
|
Activates boosting of non-nuclear neutrino losses. |
|
Boost factor for non-nuclear neutrino losses. |
|
Saves rate of change of each species due to each reaction (only for |
12. Velocity damping¶
Option |
Meaning |
|---|---|
|
Enables velocity damping. All momentum components are damped thanks to a rhs source term of the form \(-\nu f(r_\mathrm{min},r_\mathrm{max}) \rho \mathbf{u}\), where \(\nu\) ( |
|
Maximum damping frequency in \(\mathrm{s}^-1\) |
|
Radius at which damping starts. |
|
Radius beyond which the damping is set to its maximum value. |
|
Enables a time-dependent velocity damping. Only works with |
|
Time in \(\mathrm{s}\) of simulation time for which to apply the full maximum damping frequency. |
|
Time in \(\mathrm{s}\) of simulation time when velocity damping is set to |
13. Simulation stopping criteria¶
Option |
Meaning |
|---|---|
|
The maximum simulation time. |
|
The maximum number of integration steps. |
|
Maximum wall clock time of the job (in seconds). |
|
If enabled, the simulation will be interrupted if the maximum temperature on the grid exceeds the value provided in the |
|
Maximum temperature beyond which the simulation will be interrupted. |
14. Output options¶
Option |
Meaning |
|---|---|
|
Enables output dumping every |
|
Dumping frequency (only if |
|
Enables output dumping every (approximately) |
|
Dumping frequency (only if |
|
Frequency at which the simulation time, cycle, and time step are written to terminal. Further information on the gravity solver or the thermal diffusion step (if used) is also shown. |
|
Dumps a restart file every |
|
Restarts the simulation from the restart file indicated in the last line of |
|
Dumps 2D slices (planes) to a separate hdf5 output channel (only if |
|
Number of |
|
Number of |
|
Number of |
|
Output cadence of the planes in units of simulation time. |
|
Dumps spherical projections onto given radii (only if |
|
number of spherical projections |
|
Output cadence of the spherical projections in units of simulation time. In |
|
This option allows the value of state variables in certain cells to be saved to file at very timestep. The number of probes ( |
|
The number of point probes used in the run. |
|
If >1, the code performs rebinning to save a grid snapshot to the output at reduced resolution according to the provided compression factor (CF), without modifying the restart files (only for |
13. Specs¶
Option |
Meaning |
|---|---|
|
floating-point precision mode |
|
endianness setting for binary I/O paths |
|
enforces mpi-barrier calls before communication calls |