# Hydrodynamics Units

The `Hydro`

unit solves Euler’s equations for compressible gas
dynamics in one, two, or three spatial dimensions. We first describe the
basic functionality; see implementation sections below for various
extensions.

The Euler equations can be written in conservative form as

where \(\rho\) is the fluid density, \({\bf v}\) is the fluid velocity, \(P\) is the pressure, \(E\) is the sum of the internal energy \(\epsilon\) and kinetic energy per unit mass,

\({\bf g}\) is the acceleration due to gravity, and \(t\) is the time coordinate. The pressure is obtained from the energy and density using the equation of state. For the case of an ideal gas equation of state, the pressure is given by

where \(\gamma\) is the ratio of specific heats. More general equations of state are discussed in and .

In regions where the kinetic energy greatly dominates the total energy, computing the internal energy using

can lead to unphysical values, primarily due to truncation error. This results in inaccurate pressures and temperatures. To avoid this problem, we can separately evolve the internal energy according to

If the internal energy is a small fraction of the kinetic energy
(determined via the runtime parameter `Eos/eintSwitch`

), then the
total energy is recomputed using the internal energy from
[Eqn:evolve_eint] and the velocities from the
momentum equation. Numerical experiments using the PPM solver included
with Flash-X showed that using [Eqn:evolve_eint]
when the internal energy falls below \(10^{-4}\) of the kinetic
energy helps avoid the truncation errors while not affecting the
dynamics of the simulation.

For reactive flows, a separate advection equation must be solved for each chemical or nuclear species

where \(X_\ell\) is the mass fraction of the \(\ell\)th
species, with the constraint that \(\sum_\ell X_\ell = 1\). Flash-X
will enforce this constraint if you set the runtime parameter
`irenorm`

equal to 1. Otherwise, Flash-X will only restrict the
abundances to fall between `smallx`

and 1. The quantity
\(\rho X_\ell\) represents the partial density of the
\(\ell\)th fluid. The code does not explicitly track interfaces
between the fluids, so a small amount of numerical mixing can be
expected during the course of a calculation.

The `hydro`

unit has a capability to advect mass scalars. Mass scalars
are field variables advected with density, similar to species mass
fractions,

where \(\phi_\ell\) is the \(\ell\)th mass scalar. Note that
mass scalars are optional variables; to include them specify the name of
each mass scalar in a `Config`

file using the `MASS_SCALAR`

keyword.
Mass scalars are not renormalized in order to sum to 1, except when they
are declared to be part of a renormalization group. See for more
details.

## Gas hydrodynamics

### Usage

The two gas hydrodynamic solvers supplied in the release of Flash-X are organized into two different operator splitting methods: directionally split and unsplit. The directionally split piecewise-parabolic method (PPM) makes use of second-order Strang time splitting, and the new directionally unsplit solver is based on Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) Hancock type second-order scheme.

The algorithms are described in and and implemented in the directory
tree under `physics/Hydro/HydroMain/split/PPM`

and
`physics/Hydro/HydroMain/unsplit/Hydro_Unsplit`

.

Current and future implementations of `Hydro`

use the runtime
parameters and solution variables described in and . Additional runtime
parameters used either solely by the PPM method or the unsplit hydro
solver are described in `Hydro/HydroMain`

.

(`Hydro`

) unit.

Variable

Type

Default

Description

`eintSwitch`

real

0

If :math:`epsil on < {tt eintSwitch}

cdot {1

over 2}|{bf v}|^2`, use the internal energy equation to update the pressure

`irenorm`

integer

0

If equal to one, renormalize multifluid abundances following a hydro update; else restrict their values to lie between

`smallx`

and 1.

`cfl`

real

0.8

Co urant-Friedrichs-Lewy (CFL) factor; must be less than 1 for stability in explicit schemes

(`Hydro`

) unit.

Variable

Type

Description

`dens`

PER_VOLUME

density

`velx`

PER_MASS

\(x\)-component of velocity

`vely`

PER_MASS

\(y\)-component of velocity

`velz`

PER_MASS

\(z\)-component of velocity

`pres`

GENERIC

pressure

`ener`

PER_MASS

specific total energy (\(T+U\))

`temp`

GENERIC

temperature

### The unsplit hydro solver

A directionally unsplit pure hydrodynamic solver (unsplit hydro) is an alternate gas dynamics solver to the split PPM scheme. The method basically adopts a predictor-corrector type formulation (zone-edge data-extrapolated method) that provides second-order solution accuracy for smooth flows and first-order accuracy for shock flows in both space and time. Recently, the order of spatial accuracy in data reconstruction for the normal direction has been extended to implement the 3rd order PPM and 5th order Weighted ENO (WENO) methods. This unsplit hydro solver can be considered as a reduced version of the Unsplit Staggered Mesh (USM) MHD solver (see details in ) that has been available in previous Flash-X releases.

The unsplit hydro implementation can solve 1D, 2D and 3D problems with added capabilities of exploring various numerical implementations: different types of Riemann solvers; slope limiters; first, second, third and fifth reconstruction methods; a strong shock/rarefaction detection algorithm as well as two different entropy fix routines for Roe’s linearized Riemann solver.

One of the notable features of the unsplit hydro scheme is that it particularly improves the preservation of flow symmetries as compared to the splitting formulation. Also, the scheme used in this unsplit algorithm can take a wide range of CFL stability limits (e.g., CFL \(<\) 1) for all three dimensions, which is based on using upwinded transverse flux formulations developed in the multidimensional USM MHD solver (Lee, 2006; Lee and Deane, 2009; Lee, 2013).

Schemes* in the unsplit hydro solver
(`physics/Hydro/HydroMain/unsplit/Hydro_Unsplit`

)

Variable

Type

Default

Description

`order`

integer

2

Order of method in data reconstruction: 1st order Godunov (FOG), 2nd order MUSCL-Hancock (MH), 3rd order PPM, 5th order WENO.

`transOrder`

integer

1

Interpolation order of accuracy of taking upwind biased transverse flux derivatives in the unsplit data reconstruction: 1st, 2nd, 3rd. The choice of using

`transOrder`

=4 adopts a slope limiter between the 1st and 3rd order accurate methods to minimize oscillations in upwinding at discontinuities.

`slopeLimiter`

string

“vanLeer”

Slope limiter: “MINMOD”, “MC”, “VANLEER”, “HYBRID”, “LIMITED”

`LimitedSlopeBeta`

real

1.0

Slope parameter specific for the “LIMITED” slope by Toro

`charLimiting`

logical

.true.

Enable/disable limiting on characteristic variables (.false. will use limiting on primitive variables)

`use_steepening`

logical

.false.

Enable/disable contact discontinuity steepening for PPM and WENO

`use_flattening`

logical

.false.

Enable/disable flattening (or reducing) numerical oscillations for MH, PPM, and WENO

`use_avisc`

logical

.false.

Enable/disable artificial viscosity for FOG, MH, PPM, and WENO

`cvisc`

real

0.1

Artificial viscosity coefficient

`use_upwindTVD`

logical

.false.

Enable/disable upwinded TVD slope limiter PPM. NOTE: This requires NGUARD=6

`use_hybridOrder`

logical

.false.

Enable an adaptively varying reconstruction order scheme reducing its order from a high-order to first-order depending on monotonicity constraints

`` use_gravHalfUpdate``

logical

.false.

On/off gravitational acceleration source terms at the half time Riemann state update

`use_3dFullCTU`

logical

.true.

Enable a full CTU (e.g., similar to the standard 12-Riemann solve) algorithm that provides full CFL stability in 3D. If .false., then the theoretical CFL bound for 3D becomes less than 0.5 based on the linear Fourier analysis.

the unsplit hydro solver
(`physics/Hydro/HydroMain/unsplit/Hydro_Unsplit`

)

Variable

Type

Default

Description

` RiemannSolver`

string

“Roe”

Different choices for Riemann solver. “LLF (local L ax-Friedrichs)”, “HLL”, “HLLC”, “HYBRID”, “ROE”, and “Marquina”

`shockDetect`

logical

.false.

On/off attempting to detect strong sho cks/rarefactions (and saving flag in

`"shok"`

variable)` shockLowerCFL`

logical

.false.

On/off lowering of CFL factor where strong shocks are detected, automatically sets

`shockDetect`

if on.` EOSforRiemann`

logical

.false.

Enable/disable calling EOS in computing each Godunov flux

`entropy`

logical

.false.

On/off entropy fix algorithm for Roe solver

`en tropyFixMethod`

string

` “HARTENHYMAN”`

Entropy fix method for the Roe solver. “HARTEN”, “HARTENHYMAN”

The above set of runtime parameters provide various types of different combinations that help in obtaining numerical accuracy, efficiency and stability. However, there are some important tips users should know before using them.

[Extended stencil]: When

`NGUARD=6`

is used, users should also use`nxb, nyb`

, and`nzb`

larger than`2*NGUARD`

. For example, specifying`-nxb=16`

in the setup works well for 1D cases. Once setting up`NGUARD=6`

, users still can use FOG, MH, PPM, or WENO without changing`NGUARD`

back to 4.[

`transOrder`

]: The first order method`transOrder=1`

is a default and only supported method that is stable according to the linear Fourier stability analysis. The choices for higher-order interpolations are no longer available in this release.[

`EOSforRiemann`

]:`EOSforRiemann = .true.`

will call (expensive) EOS routines to compute consistent adiabatic indices (*i.e.*,`gamc, game`

) according to the given left and right states in Riemann solvers. For the ideal gamma law, in which those adiabatic indices are constant, it is not required to call EOS at all and users can set it`.false.`

to reduce computation time. On the other hand, for a degenerate gas, one can enable this switch to compute thermodynamically consistent`gamc, game`

, which in turn are used to compute the sound speed and internal energy in Riemann flux calculations. When disabled, interpolations will be used instead to get approximations of`gamc, game`

. This interpolation method has been tested and proven to gain significant computational efficiency and accuracy, giving reliable numerical solutions even for simulating a degenerate gas.[Gravity coupling with Unplit Hydro Solvers]: When gravity is included in a simulation using the unsplit hydro and MHD solvers, users can choose to include gravitational source terms in the Riemann state update at \(n+1/2\) time step (i.e.,

`use_gravHalfUpdate`

=.true.). This will provide an improved second-order accuracy with respect to coupling gravitational accelerations to hydrodynamics. It should be noted that current optimized unsplit hydro/MHD codes (*e.g.*, those selected with`+uhd, +usm`

) do not support the runtime parameters`use_gravPotUpdate`

and`use_gravConsv`

of some previous Flash-X versions any more.[Reduced CTU vs. Full CTU for 3D in the unsplit hydro (UHD) and staggered mesh (USM) solvers]:

`use_3dFullCTU`

is a new switch that enhances a numerical stability for 3D simulations in the unsplit solvers using the corner transport upwind (CTU) algorithm by Colella. The unsplit solvers of Flash-X are different from many other shock capturing codes, in that neither UHD nor USM solvers need intermediate Riemann solver solutions for updating transverse fluxes in multidimensional problems. This provides a computational efficiency because there is a reduced number of calls to Riemann solvers per cell per time step. The total number of required Riemann solver solutions are two for 2D and three for 3D (except for extra Riemann calls for constraint-transport (CT) update in USM). This is smaller than the usual stabilty requirement in many other codes which needs four for 2D and twelve for 3D in order to provide a full CFL limit (*i.e.*, CFL \(<1\)).In general for 3D, there is another computationally efficient approach that only uses six Riemann solutions (aka, 6-CTU) instead of solving twelve Riemann problems (aka, 12-CTU). In this efficient 6-CTU, however, the numerical stability limit becomes CFL\(<0.5\).

For solving 3D problems in UHD and USM, enabling the new switch

`use_3dFullCTU=.true.`

(i.e., full-CTU) will make the solution evolution scheme similar to 12-CTU while requiring to solve three Riemann problems only (again, except for the CT update in USM). On the other hand,`use_3dFullCTU=.false.`

(i.e., reduced-CTU) will be similar to the 6-CTU integration algorithm with a reduced CFL limit (*i.e.*, CFL \(<0.5\)).

The unsplit solver can take a wide range of CFL limits in all three
dimensions (*i.e.*, CFL \(<\) 1). However, in some circumstances
where there are strong shocks and rarefactions,
`shockLowerCFL=.true.`

could be useful to gain more numerical
stability by lowering the CFL accordingly (e.g., default settings
provide 0.45 for 2D and 0.25 for 3D for the Donor scheme). This
approach will automatically revert such reduced stability conditions
to any given original condition set by users when there are no
significant shocks and rarefactions detected.