compressible_flow

The dream.compressible_flow module provides a high-order HDG/DG solver for the compressible Navier-Stokes equations. It covers everything from equation-of-state and viscosity models through Riemann solvers and non-reflecting boundary conditions to a wide selection of implicit, explicit, and IMEX time integration schemes.

Note

We currently only support two-dimensional domains.

Solver Architecture

CompressibleFlowSolver is configured by composing interchangeable sub-objects for the finite element method, equation of state, viscosity model, non-dimensionalisation scaling, and Riemann solver. As with all dream solvers, solver.time selects the outer solution loop and solver.fem.scheme selects the numerical time integration scheme.

Compressible flow equations

The fundamental equations describing the motion of an unsteady, viscous and compressible flow in a space-time cylinder \(\Omega \times (0, t_{end}] \in \mathbb{R}^{d+1}\) with non-empty bounded \(d\)-dimensional spatial domain \(\Omega\), with boundary \(\partial \Omega\), and final time \(t_{end}\), are specified by the Navier-Stokes equations. In terms of conservative variables, \(\vec{U} = \begin{pmatrix} \rho, & \rho \vec{u}, & \rho E \end{pmatrix}^\T\), with density \(\rho\), velocity \(\vec{u}\), and total specific energy \(E\), this system can be expressed in dimensionless form as

\[\begin{align*} \frac{\partial \vec{U}}{\partial t} + \div(\vec{F}(\vec{U}) - \vec{G}(\vec{U}, \nabla \vec{U} )) = \vec{0}. \end{align*}\]

In a general form the convective \(\vec{F}\) and the viscous fluxes \(\vec{G}\) are given by

\[\begin{align*} \vec{F}(\vec{U}) & = \begin{pmatrix} \rho \vec{u}^\T \\ \rho \vec{u} \otimes \vec{u} + p \I \\ \rho H \vec{u}^\T \end{pmatrix}, & \vec{G}(\vec{U}, \nabla \vec{U}) = \begin{pmatrix}\vec{0}^\T \\\mat{\tau} \\ (\mat{\tau} \vec{u} - \vec{q})^\T \end{pmatrix}, \end{align*}\]

where \(p\) denotes the pressure, \(H = E + p/\rho\) the specific enthalpy, \(\mat{\tau}\) the deviatoric stress tensor, and \(\vec{q}\) the heat flux vector.

To close the system of equations we need to specify the equation of state (see eos) and the constitutive relations for the deviatoric stress tensor \(\mat{\tau}\) and the heat flux vector \(\vec{q}\) (see viscosity).

Quasi-linear Euler equations

The hyperbolic nature of the Navier-Stokes equations with respect to time lies in the Euler equations [3]

\[\begin{align*} \frac{\partial \vec{U}}{\partial t} + \div(\vec{F}(\vec{U})) = \vec{0}, \end{align*}\]

which are derived by neglecting the viscous contributions. From a characteristic point of view, it is essential to express these equations in quasi-linear form

\[\begin{align*} \frac{\partial \vec{U}}{\partial t} + \sum_{i=1}^d \mat{A}_i \frac{\partial \vec{U}}{\partial x_i} &= \vec{0}, \\ \end{align*}\]

where the \(\mat{A}_i\) are the directional convective Jacobians.

Flow Configuration

Before selecting a discretisation, the physical model is configured by setting the equation of state, viscosity model, non-dimensionalisation scaling, Mach number, and Riemann solver.

Equation of State

The equation of state closes the system by relating pressure, density, and temperature. The only currently supported option is the ideal gas law:

cfg.equation_of_state                     = 'ideal'
cfg.equation_of_state.heat_capacity_ratio = 1.4

Class	Key	Description
`IdealGas`	`'ideal'`	Ideal gas: \(p = (\gamma-1)\,\rho E_i\), \(T = \gamma E_i\)

Viscosity

The viscosity model controls the deviatoric stress tensor \(\mat{\tau}\) and heat flux \(\vec{q}\):

cfg.dynamic_viscosity = 'inviscid'   # Euler equations (no viscosity)
cfg.dynamic_viscosity = 'constant'   # constant dynamic viscosity
cfg.dynamic_viscosity = 'sutherland' # temperature-dependent viscosity

Class	Key	Description
`Inviscid`	`'inviscid'`	\(\mat{\tau} = 0\), \(\vec{q} = 0\) — Euler equations
`Constant`	`'constant'`	Constant dynamic viscosity \(\mu\)
`Sutherland`	`'sutherland'`	Sutherland’s law: \(\mu \propto T^{3/2}/(T + S)\)

Non-Dimensionalisation Scaling

The scaling determines how the reference state \((\rho_\infty, c_\infty, p_\infty, T_\infty)\) is constructed from the Mach number. It affects all non-dimensional quantities used by the solver.

cfg.scaling     = 'acoustic'      # acoustic: rho_inf=1, c_inf=1, p_inf=1/gamma
cfg.scaling     = 'aerodynamic'   # aerodynamic: rho_inf=1, u_inf=1, p_inf=1/(gamma*M^2)
cfg.scaling     = 'aeroacoustic'  # combined aerodynamic/acoustic scaling
cfg.mach_number = 0.03

Class	Key	Characteristic velocity
`Acoustic`	`'acoustic'`	\(u_\mathrm{ref} = c_\infty\)
`Aerodynamic`	`'aerodynamic'`	\(u_\mathrm{ref} = u_\infty\)
`Aeroacoustic`	`'aeroacoustic'`	\(u_\mathrm{ref} = \sqrt{u_\infty c_\infty}\)

Riemann Solver

The Riemann solver defines the convective numerical flux at element interfaces and boundary faces. It affects both accuracy and stability:

cfg.riemann_solver = 'lax_friedrich'

Class	Key	Notes
`LaxFriedrich`	`'lax_friedrich'`	Simple, robust; adds maximum-speed dissipation
`Upwind`	`'upwind'`	Characteristic upwinding
`Roe`	`'roe'`	Roe-averaged flux; accurate but no entropy fix
`HLL`	`'hll'`	Two-wave HLL flux
`HLLEM`	`'hllem'`	HLL with entropy-fix correction

Discretisation

The solver discretises the compressible Navier-Stokes equations on a mesh \(\mesh\) using either the Hybridised Discontinuous Galerkin (HDG) method (ConservativeHDG) or the standard Discontinuous Galerkin (DG) method (ConservativeDG).

HDG introduces an additional skeleton unknown \(\widehat{\vec{U}}_h\) living on the mesh facets \(\facets\). The local element problems are solved independently for each element \(K \in \mesh\) given \(\widehat{\vec{U}}_h\), and then a global system is assembled solely in terms of the skeleton degrees of freedom. This static condensation makes HDG particularly efficient for high-order discretisations since the global system is much smaller than the full DG system. The viscous terms are handled by the viscous_treatment sub-object of ConservativeHDG: the default StrainHeat formulation uses the physical strain-rate and heat-flux form; Gradient and InteriorPenaltyHDG are available as alternatives.

DG keeps all degrees of freedom element-local with no additional skeleton unknowns. It is straightforward to use with explicit time integrators, at the cost of a larger globally coupled system for implicit schemes. The viscous treatment for DG is InteriorPenaltySDG.

In both formulations the convective numerical flux at element interfaces is determined by the selected RiemannSolver. The available options range from the simple LaxFriedrich solver to the more accurate Upwind, Roe, HLL, and HLLEM solvers.

Boundary Conditions

Boundary conditions are assigned to named mesh boundaries via the solver’s bcs (BoundaryConditions) attribute. The following condition types are available:

Class	Description
`FarField`	Characteristic inflow/outflow; farfield state \((\rho, \mathbf{u}, T)_\infty\)
`Outflow`	Pressure outflow: prescribe static pressure \(p\)
`Inflow`	Inflow: prescribe \((\rho, \mathbf{u})\) or total conditions
`InviscidWall`	Slip wall: zero normal velocity
`Symmetry`	Symmetry plane
`IsothermalWall`	No-slip wall with fixed temperature
`AdiabaticWall`	No-slip wall with zero heat flux
`CBC`	Non-reflecting BC; subclasses: `GRCBC`, `NSCBC`
`InterfaceBC`	IMEX coupling: solution fields from neighbouring submesh at \(\Gamma_i\)

Initial and Domain Conditions

Domain conditions are assigned via the solver’s dcs (DomainConditions) attribute and apply over named mesh regions rather than boundaries. For transient simulations the initial state is mandatory:

Uinf = solver.get_farfield_fields((1, 0))
solver.dcs['domain'] = Initial(fields=Uinf)

The fields object passed to Initial must be a flowfields instance, typically obtained via get_farfield_fields() for a uniform free-stream state, or assembled manually from density, velocity, and pressure (or temperature) using flowfields(rho=..., u=..., p=...).

Class	Description
`Initial`	Initial state \(\mathbf{U}(\cdot,0)\) projected onto the solution space; required for `'transient'` and `'pseudo_time_stepping'` outer loops
`Force`	Body force term added to the momentum equation
`Perturbation`	Superimposed perturbation field on the background state
`SpongeLayer`	Volumetric damping term that drives the solution toward a `target_state` with a prescribed weight `function`; used to absorb outgoing waves at domain boundaries
`PSpongeLayer`	Polynomial sponge layer: same damping as `SpongeLayer` plus gradual reduction of the local polynomial order from `high_order` to `low_order` for additional dissipation
`GridDeformation`	Mesh deformation mapping defined by `GridMapping` objects per coordinate direction; deforms the mesh in a domain region

Time Integration

The outer solution loop is selected via solver.time and the numerical time integration scheme via solver.fem.scheme (see the time module for the full time infrastructure).

Implicit schemes (unconditionally stable — recommended for viscous or diffusion-dominated problems and large time steps):

Scheme	Stages	Order	Key
`ImplicitEuler`	1	1	`'implicit_euler'`
`BDF2`	1	2	`'bdf2'`
`BDF3`	1	3	`'bdf3'`
`SDIRK22`	2	2	`'sdirk22'`
`SDIRK33`	3	3	`'sdirk33'`
`SDIRK43`	4	3	`'sdirk43'`
`SDIRK54`	5	4	`'sdirk54'`
`DIRK43_WSO2`	4	3 (WSO 3)	`'dirk43_wso2'`
`DIRK34_LDD`	3	4	`'dirk34_ldd'`

Explicit schemes (require CFL condition — suitable for inviscid or convection-dominated problems):

Scheme	Stages	Order	Key
`ExplicitEuler`	1	1	`'explicit_euler'`
`RK_ARS22`	2	2	`'rk_ars22'`
`RK_ARS232`	2	2	`'rk_ars232'`
`RK_ARS33`	3	3	`'rk_ars33'`
`RK_ARS43`	4	3	`'rk_ars43'`
`SSPRK3`	3	3	`'ssprk3'`
`CRK4`	4	4	`'crk4'`

IMEX schemes (geometry splitting — two submeshes \(\mesh^{im} \cup \mesh^{ex}\)):

The implicit HDG solver runs on \(\mesh^{im}\) and the explicit DG solver on \(\mesh^{ex}\). Each solver uses its own scheme, but the stage times must be synchronised: \(c_i^{im} = \bar{c}_{i+1}^{ex}\). The ARS-family schemes are designed for exactly this pairing:

Order	Implicit scheme	Explicit scheme
2	`'sdirk22'`	`'rk_ars22'`
2	`'sdirk22'`	`'rk_ars232'`
3	`'sdirk33'`	`'rk_ars33'`
3	`'sdirk43'`	`'rk_ars43'`

API Reference

`solver`	Compressible flow solver configuration
`conservative`	Definitions of conservative methods
`eos`	Definitions of equation of states for compressible flow.
`riemann_solver`	Definitions of Riemann solvers for compressible flow.
`viscosity`	Definitions of viscous constitutive relations for compressible flow
`scaling`	Definitions of dimensionless compressible flow equations.
`config`	Definitions of boundary/domain conditions for compressible flow