Compressible Navier-Stokes mini-app

This example is located in the subdirectory examples/fluids. It solves the time-dependent Navier-Stokes equations of compressible gas dynamics in a static Eulerian three-dimensional frame using unstructured high-order finite/spectral element spatial discretizations and explicit or implicit high-order time-stepping (available in PETSc). Moreover, the Navier-Stokes example has been developed using PETSc, so that the pointwise physics (defined at quadrature points) is separated from the parallelization and meshing concerns.

Running the mini-app

The Navier-Stokes mini-app is controlled via command-line options. The following options are common among all problem types:

Table 1 Common Runtime Options

Option

Description

Default value

-ceed

CEED resource specifier

/cpu/self/opt/blocked

-test

Run in test mode

false

-compare_final_state_atol

Test absolute tolerance

1E-11

-compare_final_state_filename

Test filename

-problem

Problem to solve (advection, advection2d, density_current, or euler_vortex)

density_current

-implicit

Use implicit time integartor formulation

-degree

Polynomial degree of tensor product basis (must be >= 1)

1

-q_extra

Number of extra quadrature points

2

-viz_refine

Use regular refinement for visualization

0

-output_freq

Frequency of output, in number of steps

10

-continue

Continue from previous solution

0

-output_dir

Output directory

.

-bc_wall

Use wall boundary conditions on this list of faces

-wall_comps

An array of constrained component numbers for wall BCs

-bc_slip_x

Use slip boundary conditions, for the x component, on this list of faces

-bc_slip_y

Use slip boundary conditions, for the y component, on this list of faces

-bc_slip_z

Use slip boundary conditions, for the z component, on this list of faces

-bc_inflow

Use inflow boundary conditions on this list of faces

-bc_outflow

Use outflow boundary conditions on this list of faces

-snes_view

View PETSc SNES nonlinear solver configuration

-log_view

View PETSc performance log

-help

View comprehensive information about run-time options

For the case of a square/cubic mesh, the list of face indices to be used with -bc_wall, bc_inflow, bc_outflow and/or -bc_slip_x, -bc_slip_y, and -bc_slip_z are:

Table 2 2D Face ID Labels

PETSc Face Name

Cartesian direction

Face ID

faceMarkerBottom

-z

1

faceMarkerRight

+x

2

faceMarkerTop

+z

3

faceMarkerLeft

-x

4

Table 3 2D Face ID Labels

PETSc Face Name

Cartesian direction

Face ID

faceMarkerBottom

-z

1

faceMarkerTop

+z

2

faceMarkerFront

-y

3

faceMarkerBack

+y

4

faceMarkerRight

+x

5

faceMarkerLeft

-x

6

Advection

For testing purposes, there is a reduced mode for pure advection, which holds density \(\rho\) and momentum density \(\rho \bm u\) constant while advecting “total energy density” \(E\). These are available in 2D and 3D.

2D advection

For the 2D advection problem, the following additional command-line options are available:

Table 4 Advection2D Runtime Options

Option

Description

Default value

Unit

-rc

Characteristic radius of thermal bubble

1000

m

-units_meter

1 meter in scaled length units

1E-2

-units_second

1 second in scaled time units

1E-2

-units_kilogram

1 kilogram in scaled mass units

1E-6

-strong_form

Strong (1) or weak/integrated by parts (0) residual

0

-stab

Stabilization method (none, su, or supg)

none

-CtauS

Scale coefficient for stabilization tau (nondimensional)

0

-wind_type

Wind type in Advection (rotation or translation)

rotation

-wind_translation

Constant wind vector when -wind_type translation

1,0,0

-E_wind

Total energy of inflow wind when -wind_type translation

1E6

J

An example of the rotation mode can be run with:

./navierstokes -problem advection2d -dm_plex_box_faces 20,20 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1000,1000 -bc_wall 1,2,3,4 -wall_comps 4 -wind_type rotation -implicit -stab supg

and the translation mode with:

./navierstokes -problem advection2d -dm_plex_box_faces 20,20 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1000,1000 -units_meter 1e-4 -wind_type translation -wind_translation 1,-.5 -bc_inflow 1,2,3,4

Note the lengths in -dm_plex_box_upper are given in meters, and will be nondimensionalized according to -units_meter.

3D advection

For the 3D advection problem, the following additional command-line options are available:

Table 5 Advection3D Runtime Options

Option

Description

Default value

Unit

-rc

Characteristic radius of thermal bubble

1000

m

-units_meter

1 meter in scaled length units

1E-2

-units_second

1 second in scaled time units

1E-2

-units_kilogram

1 kilogram in scaled mass units

1E-6

-strong_form

Strong (1) or weak/integrated by parts (0) residual

0

-stab

Stabilization method (none, su, or supg)

none

-CtauS

Scale coefficient for stabilization tau (nondimensional)

0

-wind_type

Wind type in Advection (rotation or translation)

rotation

-wind_translation

Constant wind vector when -wind_type translation

1,0,0

-E_wind

Total energy of inflow wind when -wind_type translation

1E6

J

-bubble_type

sphere (3D) or cylinder (2D)

shpere

-bubble_continuity

smooth, back_sharp, or thick

smooth

An example of the rotation mode can be run with:

./navierstokes -problem advection -dm_plex_box_faces 10,10,10 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 8000,8000,8000 -bc_wall 1,2,3,4,5,6 -wall_comps 4 -wind_type rotation -implicit -stab su

and the translation mode with:

./navierstokes -problem advection -dm_plex_box_faces 10,10,10 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 8000,8000,8000 -wind_type translation -wind_translation .5,-1,0 -bc_inflow 1,2,3,4,5,6

Inviscid Ideal Gas

Isentropic Euler vortex

For the Isentropic Vortex problem, the following additional command-line options are available:

Table 6 Isentropic Vortex Runtime Options

Option

Description

Default value

Unit

-center

Location of vortex center

(lx,ly,lz)/2

(m,m,m)

-units_meter

1 meter in scaled length units

1E-2

-units_second

1 second in scaled time units

1E-2

-mean_velocity

Background velocity vector

(1,1,0)

-vortex_strength

Strength of vortex < 10

5

-c_tau

Stabilization constant

0.5

This problem can be run with:

./navierstokes -problem euler_vortex -dm_plex_box_faces 20,20,1 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1000,1000,50 -dm_plex_dim 3 -bc_inflow 4,6 -bc_outflow 3,5 -bc_slip_z 1,2 -mean_velocity .5,-.8,0.

Sod shock tube

For the Shock Tube problem, the following additional command-line options are available:

Table 7 Shock Tube Runtime Options

Option

Description

Default value

Unit

-units_meter

1 meter in scaled length units

1E-2

-units_second

1 second in scaled time units

1E-2

-yzb

Use YZB discontinuity capturing

none

-stab

Stabilization method (none, su, or supg)

none

This problem can be run with:

./navierstokes -problem shocktube -yzb -stab su -bc_slip_z 3,4 -bc_slip_y 1,2 -bc_wall 5,6 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1000,100,100 -dm_plex_box_faces 200,1,1 -units_second 0.1 

Newtonian viscosity, Ideal Gas

For the Density Current, Channel, and Blasius problems, the following common command-line options are available:

Table 8 Newtonian Ideal Gas problems Runtime Options

Option

Description

Default value

Unit

-units_meter

1 meter in scaled length units

1

-units_second

1 second in scaled time units

1

-units_kilogram

1 kilogram in scaled mass units

1

-units_Kelvin

1 Kelvin in scaled temperature units

1

-stab

Stabilization method (none, su, or supg)

none

-c_tau

Stabilization constant, \(c_\tau\)

0.5

-Ctau_t

Stabilization time constant, \(C_t\)

1.0

-Ctau_v

Stabilization viscous constant, \(C_v\)

36.0

-Ctau_C

Stabilization continuity constant, \(C_c\)

1.0

-Ctau_M

Stabilization momentum constant, \(C_m\)

1.0

-Ctau_E

Stabilization energy constant, \(C_E\)

1.0

-cv

Heat capacity at constant volume

717

J/(kg K)

-cp

Heat capacity at constant pressure

1004

J/(kg K)

-g

Gravitational acceleration

9.81

m/s^2

-lambda

Stokes hypothesis second viscosity coefficient

-2/3

-mu

Shear dynamic viscosity coefficient

75

Pa s

-k

Thermal conductivity

0.02638

W/(m K)

-newtonian_unit_tests

Developer option to test properties

false

boolean

Density current

The Density Current problem the following command-line options are available in addition to the Newtonian Ideal Gas options:

Table 9 Density Current Runtime Options

Option

Description

Default value

Unit

-center

Location of bubble center

(lx,ly,lz)/2

(m,m,m)

-dc_axis

Axis of density current cylindrical anomaly, or (0,0,0) for spherically symmetric

(0,0,0)

-rc

Characteristic radius of thermal bubble

1000

m

-theta0

Reference potential temperature

300

K

-thetaC

Perturbation of potential temperature

-15

K

-P0

Atmospheric pressure

1E5

Pa

-N

Brunt-Vaisala frequency

0.01

1/s

This problem can be run with:

./navierstokes -problem density_current -dm_plex_box_faces 16,1,8 -degree 1 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 2000,125,1000 -dm_plex_dim 3 -rc 400. -bc_wall 1,2,5,6 -wall_comps 1,2,3 -bc_slip_y 3,4 -mu 75

Channel flow

The Channel problem the following command-line options are available in addition to the Newtonian Ideal Gas options:

Table 10 Channel Runtime Options

Option

Description

Default value

Unit

-umax

Maximum/centerline velocity of the flow

10

m/s

-theta0

Reference potential temperature

300

K

-P0

Atmospheric pressure

1E5

Pa

-body_force_scale

Multiplier for body force (-1 for flow reversal)

1

This problem can be run with the channel.yaml file via:

./navierstokes -options_file channel.yaml
problem: 'channel'
mu: .01

umax: 40
implicit: true
ts:
  type: 'beuler'
  adapt_type: 'none'
  dt: 5e-6

dm_plex_box_lower: 0,0,0
dm_plex_box_upper: .01,.01,.001
dm_plex_dim: 3
degree: 1
dm_plex_box_faces: 10,10,1
bc_slip_z: 1,2
bc_wall: 3,4
wall_comps: 1,2,3
dm_plex_box_bd: 'periodic,none,none'

Blasius boundary layer

The Blasius problem the following command-line options are available in addition to the Newtonian Ideal Gas options:

Table 11 Blasius Runtime Options

Option

Description

Default value

Unit

-Uinf

Freestream velocity

40

m/s

-delta0

Boundary layer height at the inflow

4.2e-4

m

-theta0

Reference potential temperature

288

K

-P0

Atmospheric pressure

1.01E5

Pa

-platemesh_refine_height

Height at which -platemesh_Ndelta number of elements should refined into

5.9E-4

m

-platemesh_Ndelta

Number of elements to keep below -platemesh_refine_height

45

-platemesh_growth

Growth rate of the elements in the refinement region

1.08

-platemesh_top_angle

Downward angle of the top face of the domain. This face serves as an outlet.

5

degrees

-stg_use

Whether to use stg for the inflow conditions

false

-platemesh_y_node_locs_path

Path to file with y node locations. If empty, will use mesh warping instead.

""

This problem can be run with the blasius.yaml file via:

./navierstokes -options_file blasius.yaml
problem: 'blasius'

implicit: true
ts:
  adapt_type: 'none'
  type: 'beuler'
  dt: 0.2e-5
  max_time: 1.0e-3
output_freq: 10

## Linear Settings:
degree: 1
dm_plex_box_faces: 40,60,1
platemesh_nDelta: 45

# # Quadratic Settings:
# degree: 2
# dm_plex_box_faces: 20,30,1
# platemesh:
#   nDelta: 22
#   growth: 1.1664 # 1.08^2

stab: 'supg'
Ctau_t: 1
#Ctau_v: 36,60,128 is what PHASTA has for p=1,2,3
# Linear Settings:
Ctau_v: 36
Ctau_C: 0.25
Ctau_M: 0.25
Ctau_E: 0.125
# # Quadratic Settings:
# Ctau_v: 60
# Ctau_C: 0.125
# Ctau_M: 0.125
# Ctau_E: 0.125

q_extra: 0

dm_plex_box_lower: 0,0,0
dm_plex_box_upper: 4.2e-3,4.2e-3,5.e-5
dm_plex_dim: 3
# Faces labeled 1=z- 2=z+ 3=y- 4=y+ 5=x+ 6=x-
bc_slip_z: 1,2
bc_wall: 3
wall_comps: 1,2,3
bc_inflow: 6
bc_outflow: 5,4
g: 0,0,0

stg:
  use: false
  inflow_path: "./STGInflow_blasius.dat"
  mean_only: true

STG Inflow for Flat Plate

Using the STG Inflow for the blasius problem adds the following command-line options:

Table 12 Blasius Runtime Options

Option

Description

Default value

Unit

-stg_inflow_path

Path to the STGInflow file

./STGInflow.dat

-stg_rand_path

Path to the STGRand file

./STGRand.dat

-stg_alpha

Growth rate of the wavemodes

1.01

-stg_u0

Convective velocity, \(U_0\)

0.0

m/s

-stg_mean_only

Only impose the mean velocity (no fluctutations)

false

-stg_strong

Strongly enforce the STG inflow boundary condition

false

This problem can be run with the blasius.yaml file via:

./navierstokes -options_file blasius.yaml -stg_use true

Note the added -stg_use true flag This overrides the stg: use: false setting in the blasius.yaml file, enabling the use of the STG inflow.

The Navier-Stokes equations

The mathematical formulation (from [GRL10], cf. SE3) is given in what follows. The compressible Navier-Stokes equations in conservative form are

(15)\[ \begin{aligned} \frac{\partial \rho}{\partial t} + \nabla \cdot \bm{U} &= 0 \\ \frac{\partial \bm{U}}{\partial t} + \nabla \cdot \left( \frac{\bm{U} \otimes \bm{U}}{\rho} + P \bm{I}_3 -\bm\sigma \right) + \rho g \bm{\hat k} &= 0 \\ \frac{\partial E}{\partial t} + \nabla \cdot \left( \frac{(E + P)\bm{U}}{\rho} -\bm{u} \cdot \bm{\sigma} - k \nabla T \right) &= 0 \, , \\ \end{aligned} \]

where \(\bm{\sigma} = \mu(\nabla \bm{u} + (\nabla \bm{u})^T + \lambda (\nabla \cdot \bm{u})\bm{I}_3)\) is the Cauchy (symmetric) stress tensor, with \(\mu\) the dynamic viscosity coefficient, and \(\lambda = - 2/3\) the Stokes hypothesis constant. In equations (15), \(\rho\) represents the volume mass density, \(U\) the momentum density (defined as \(\bm{U}=\rho \bm{u}\), where \(\bm{u}\) is the vector velocity field), \(E\) the total energy density (defined as \(E = \rho e\), where \(e\) is the total energy), \(\bm{I}_3\) represents the \(3 \times 3\) identity matrix, \(g\) the gravitational acceleration constant, \(\bm{\hat{k}}\) the unit vector in the \(z\) direction, \(k\) the thermal conductivity constant, \(T\) represents the temperature, and \(P\) the pressure, given by the following equation of state

(16)\[ P = \left( {c_p}/{c_v} -1\right) \left( E - {\bm{U}\cdot\bm{U}}/{(2 \rho)} - \rho g z \right) \, , \]

where \(c_p\) is the specific heat at constant pressure and \(c_v\) is the specific heat at constant volume (that define \(\gamma = c_p / c_v\), the specific heat ratio).

The system (15) can be rewritten in vector form

(17)\[ \frac{\partial \bm{q}}{\partial t} + \nabla \cdot \bm{F}(\bm{q}) -S(\bm{q}) = 0 \, , \]

for the state variables 5-dimensional vector

\[ \bm{q} = \begin{pmatrix} \rho \\ \bm{U} \equiv \rho \bm{ u }\\ E \equiv \rho e \end{pmatrix} \begin{array}{l} \leftarrow\textrm{ volume mass density}\\ \leftarrow\textrm{ momentum density}\\ \leftarrow\textrm{ energy density} \end{array} \]

where the flux and the source terms, respectively, are given by

(18)\[ \begin{aligned} \bm{F}(\bm{q}) &= \underbrace{\begin{pmatrix} \bm{U}\\ {(\bm{U} \otimes \bm{U})}/{\rho} + P \bm{I}_3 \\ {(E + P)\bm{U}}/{\rho} \end{pmatrix}}_{\bm F_{\text{adv}}} + \underbrace{\begin{pmatrix} 0 \\ - \bm{\sigma} \\ - \bm{u} \cdot \bm{\sigma} - k \nabla T \end{pmatrix}}_{\bm F_{\text{diff}}},\\ S(\bm{q}) &= - \begin{pmatrix} 0\\ \rho g \bm{\hat{k}}\\ 0 \end{pmatrix}. \end{aligned} \]

Let the discrete solution be

\[ \bm{q}_N (\bm{x},t)^{(e)} = \sum_{k=1}^{P}\psi_k (\bm{x})\bm{q}_k^{(e)} \]

with \(P=p+1\) the number of nodes in the element \(e\). We use tensor-product bases \(\psi_{kji} = h_i(X_0)h_j(X_1)h_k(X_2)\).

For the time discretization, we use two types of time stepping schemes.

  • Explicit time-stepping method

    The following explicit formulation is solved with the adaptive Runge-Kutta-Fehlberg (RKF4-5) method by default (any explicit time-stepping scheme available in PETSc can be chosen at runtime)

    \[ \bm{q}_N^{n+1} = \bm{q}_N^n + \Delta t \sum_{i=1}^{s} b_i k_i \, , \]

    where

    \[ \begin{aligned} k_1 &= f(t^n, \bm{q}_N^n)\\ k_2 &= f(t^n + c_2 \Delta t, \bm{q}_N^n + \Delta t (a_{21} k_1))\\ k_3 &= f(t^n + c_3 \Delta t, \bm{q}_N^n + \Delta t (a_{31} k_1 + a_{32} k_2))\\ \vdots&\\ k_i &= f\left(t^n + c_i \Delta t, \bm{q}_N^n + \Delta t \sum_{j=1}^s a_{ij} k_j \right)\\ \end{aligned} \]

    and with

    \[ f(t^n, \bm{q}_N^n) = - [\nabla \cdot \bm{F}(\bm{q}_N)]^n + [S(\bm{q}_N)]^n \, . \]
  • Implicit time-stepping method

    This time stepping method which can be selected using the option -implicit is solved with Backward Differentiation Formula (BDF) method by default (similarly, any implicit time-stepping scheme available in PETSc can be chosen at runtime). The implicit formulation solves nonlinear systems for \(\bm q_N\):

    (19)\[ \bm f(\bm q_N) \equiv \bm g(t^{n+1}, \bm{q}_N, \bm{\dot{q}}_N) = 0 \, , \]

    where the time derivative \(\bm{\dot q}_N\) is defined by

    \[ \bm{\dot{q}}_N(\bm q_N) = \alpha \bm q_N + \bm z_N \]

    in terms of \(\bm z_N\) from prior state and \(\alpha > 0\), both of which depend on the specific time integration scheme (backward difference formulas, generalized alpha, implicit Runge-Kutta, etc.). Each nonlinear system (19) will correspond to a weak form, as explained below. In determining how difficult a given problem is to solve, we consider the Jacobian of (19),

    \[ \frac{\partial \bm f}{\partial \bm q_N} = \frac{\partial \bm g}{\partial \bm q_N} + \alpha \frac{\partial \bm g}{\partial \bm{\dot q}_N}. \]

    The scalar “shift” \(\alpha\) scales inversely with the time step \(\Delta t\), so small time steps result in the Jacobian being dominated by the second term, which is a sort of “mass matrix”, and typically well-conditioned independent of grid resolution with a simple preconditioner (such as Jacobi). In contrast, the first term dominates for large time steps, with a condition number that grows with the diameter of the domain and polynomial degree of the approximation space. Both terms are significant for time-accurate simulation and the setup costs of strong preconditioners must be balanced with the convergence rate of Krylov methods using weak preconditioners.

To obtain a finite element discretization, we first multiply the strong form (17) by a test function \(\bm v \in H^1(\Omega)\) and integrate,

\[ \int_{\Omega} \bm v \cdot \left(\frac{\partial \bm{q}_N}{\partial t} + \nabla \cdot \bm{F}(\bm{q}_N) - \bm{S}(\bm{q}_N) \right) \,dV = 0 \, , \; \forall \bm v \in \mathcal{V}_p\,, \]

with \(\mathcal{V}_p = \{ \bm v(\bm x) \in H^{1}(\Omega_e) \,|\, \bm v(\bm x_e(\bm X)) \in P_p(\bm{I}), e=1,\ldots,N_e \}\) a mapped space of polynomials containing at least polynomials of degree \(p\) (with or without the higher mixed terms that appear in tensor product spaces).

Integrating by parts on the divergence term, we arrive at the weak form,

(20)\[ \begin{aligned} \int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) \,dV - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\ + \int_{\partial \Omega} \bm v \cdot \bm{F}(\bm q_N) \cdot \widehat{\bm{n}} \,dS &= 0 \, , \; \forall \bm v \in \mathcal{V}_p \,, \end{aligned} \]

where \(\bm{F}(\bm q_N) \cdot \widehat{\bm{n}}\) is typically replaced with a boundary condition.

Note

The notation \(\nabla \bm v \!:\! \bm F\) represents contraction over both fields and spatial dimensions while a single dot represents contraction in just one, which should be clear from context, e.g., \(\bm v \cdot \bm S\) contracts over fields while \(\bm F \cdot \widehat{\bm n}\) contracts over spatial dimensions.

We solve (20) using a Galerkin discretization (default) or a stabilized method, as is necessary for most real-world flows.

Galerkin methods produce oscillations for transport-dominated problems (any time the cell Péclet number is larger than 1), and those tend to blow up for nonlinear problems such as the Euler equations and (low-viscosity/poorly resolved) Navier-Stokes, in which case stabilization is necessary. Our formulation follows [HST10], which offers a comprehensive review of stabilization and shock-capturing methods for continuous finite element discretization of compressible flows.

  • SUPG (streamline-upwind/Petrov-Galerkin)

    In this method, the weighted residual of the strong form (17) is added to the Galerkin formulation (20). The weak form for this method is given as

    (21)\[ \begin{aligned} \int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) \,dV - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\ + \int_{\partial \Omega} \bm v \cdot \bm{F}(\bm{q}_N) \cdot \widehat{\bm{n}} \,dS & \\ + \int_{\Omega} \mathcal{P}(\bm v)^T \, \left( \frac{\partial \bm{q}_N}{\partial t} \, + \, \nabla \cdot \bm{F} \, (\bm{q}_N) - \bm{S}(\bm{q}_N) \right) \,dV &= 0 \, , \; \forall \bm v \in \mathcal{V}_p \end{aligned} \]

    This stabilization technique can be selected using the option -stab supg.

  • SU (streamline-upwind)

    This method is a simplified version of SUPG (21) which is developed for debugging/comparison purposes. The weak form for this method is

    (22)\[ \begin{aligned} \int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) \,dV - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\ + \int_{\partial \Omega} \bm v \cdot \bm{F}(\bm{q}_N) \cdot \widehat{\bm{n}} \,dS & \\ + \int_{\Omega} \mathcal{P}(\bm v)^T \, \nabla \cdot \bm{F} \, (\bm{q}_N) \,dV & = 0 \, , \; \forall \bm v \in \mathcal{V}_p \end{aligned} \]

    This stabilization technique can be selected using the option -stab su.

In both (22) and (21), \(\mathcal P\) is called the perturbation to the test-function space, since it modifies the original Galerkin method into SUPG or SU schemes. It is defined as

(23)\[ \mathcal P(\bm v) \equiv \bm{\tau} \left(\frac{\partial \bm{F}_{\text{adv}} (\bm{q}_N)}{\partial \bm{q}_N} \right) \, \nabla \bm v\,, \]

where parameter \(\bm{\tau} \in \mathbb R^{3}\) (spatial index) or \(\bm \tau \in \mathbb R^{5\times 5}\) (field indices) is an intrinsic time scale matrix. Most generally, we consider \(\bm\tau \in \mathbb R^{3,5,5}\). This expression contains the advective flux Jacobian, which may be thought of as mapping from a 5-vector (state) to a \((5,3)\) tensor (flux) or from a \((5,3)\) tensor (gradient of state) to a 5-vector (time derivative of state); the latter is used in (23) because it’s applied to \(\nabla\bm v\). The forward variational form can be readily expressed by differentiating \(\bm F_{\text{adv}}\) of (18)

\[ \begin{aligned} \diff\bm F_{\text{adv}}(\diff\bm q; \bm q) &= \frac{\partial \bm F_{\text{adv}}}{\partial \bm q} \diff\bm q \\ &= \begin{pmatrix} \diff\bm U \\ (\diff\bm U \otimes \bm U + \bm U \otimes \diff\bm U)/\rho - (\bm U \otimes \bm U)/\rho^2 \diff\rho + \diff P \bm I_3 \\ (E + P)\diff\bm U/\rho + (\diff E + \diff P)\bm U/\rho - (E + P) \bm U/\rho^2 \diff\rho \end{pmatrix}, \end{aligned} \]

where \(\diff P\) is defined by differentiating (16). This action is also readily computed by forward-mode AD, but since \(\bm v\) is a test function, we actually need the action of the adjoint to use (23) in finite element computation; that can be computed by reverse-mode AD. We may equivalently write the stabilization term as

\[ \mathcal P(\bm v)^T \bm r = \nabla \bm v \tcolon \left(\frac{\partial \bm F_{\text{adv}}}{\partial \bm q}\right)^T \, \bm\tau \bm r, \]

where \(\bm r\) is the strong form residual and \(\bm\tau\) is a \(5\times 5\) matrix.

Currently, this demo provides three types of problems/physical models that can be selected at run time via the option -problem. Advection, the problem of the transport of energy in a uniform vector velocity field, Isentropic Vortex, the exact solution to the Euler equations, and the so called Density Current problem.

Advection

A simplified version of system (15), only accounting for the transport of total energy, is given by

(30)\[ \frac{\partial E}{\partial t} + \nabla \cdot (\bm{u} E ) = 0 \, , \]

with \(\bm{u}\) the vector velocity field. In this particular test case, a blob of total energy (defined by a characteristic radius \(r_c\)) is transported by two different wind types.

  • Rotation

    In this case, a uniform circular velocity field transports the blob of total energy. We have solved (30) applying zero energy density \(E\), and no-flux for \(\bm{u}\) on the boundaries.

  • Translation

    In this case, a background wind with a constant rectilinear velocity field, enters the domain and transports the blob of total energy out of the domain.

    For the inflow boundary conditions, a prescribed \(E_{wind}\) is applied weakly on the inflow boundaries such that the weak form boundary integral in (20) is defined as

    \[ \int_{\partial \Omega_{inflow}} \bm v \cdot \bm{F}(\bm q_N) \cdot \widehat{\bm{n}} \,dS = \int_{\partial \Omega_{inflow}} \bm v \, E_{wind} \, \bm u \cdot \widehat{\bm{n}} \,dS \, , \]

    For the outflow boundary conditions, we have used the current values of \(E\), following [PMK92] which extends the validity of the weak form of the governing equations to the outflow instead of replacing them with unknown essential or natural boundary conditions. The weak form boundary integral in (20) for outflow boundary conditions is defined as

    \[ \int_{\partial \Omega_{outflow}} \bm v \cdot \bm{F}(\bm q_N) \cdot \widehat{\bm{n}} \,dS = \int_{\partial \Omega_{outflow}} \bm v \, E \, \bm u \cdot \widehat{\bm{n}} \,dS \, , \]

Isentropic Vortex

Three-dimensional Euler equations, which are simplified and nondimensionalized version of system (15) and account only for the convective fluxes, are given by

(31)\[ \begin{aligned} \frac{\partial \rho}{\partial t} + \nabla \cdot \bm{U} &= 0 \\ \frac{\partial \bm{U}}{\partial t} + \nabla \cdot \left( \frac{\bm{U} \otimes \bm{U}}{\rho} + P \bm{I}_3 \right) &= 0 \\ \frac{\partial E}{\partial t} + \nabla \cdot \left( \frac{(E + P)\bm{U}}{\rho} \right) &= 0 \, , \\ \end{aligned} \]

Following the setup given in [ZZS11], the mean flow for this problem is \(\rho=1\), \(P=1\), \(T=P/\rho= 1\) (Specific Gas Constant, \(R\), is 1), and \(\bm{u}=(u_1,u_2,0)\) while the perturbation \(\delta \bm{u}\), and \(\delta T\) are defined as

\[ \begin{aligned} (\delta u_1, \, \delta u_2) &= \frac{\epsilon}{2 \pi} \, e^{0.5(1-r^2)} \, (-\bar{y}, \, \bar{x}) \, , \\ \delta T &= - \frac{(\gamma-1) \, \epsilon^2}{8 \, \gamma \, \pi^2} \, e^{1-r^2} \, , \\ \end{aligned} \]

where \((\bar{x}, \, \bar{y}) = (x-x_c, \, y-y_c)\), \((x_c, \, y_c)\) represents the center of the domain, \(r^2=\bar{x}^2 + \bar{y}^2\), and \(\epsilon\) is the vortex strength (\(\epsilon\) < 10). There is no perturbation in the entropy \(S=P/\rho^\gamma\) (\(\delta S=0)\).

Shock Tube

This test problem is based on Sod’s Shock Tube (from[sod]), a canonical test case for discontinuity capturing in one dimension. For this problem, the three-dimensional Euler equations are formulated exactly as in the Isentropic Vortex problem. The default initial conditions are \(P=1\), \(\rho=1\) for the driver section and \(P=0.1\), \(\rho=0.125\) for the driven section. The initial velocity is zero in both sections. Slip boundary conditions are applied to the side walls and wall boundary conditions are applied at the end walls.

SU upwinding and discontinuity capturing have been implemented into the explicit timestepping operator for this problem. Discontinuity capturing is accomplished using a modified version of the \(YZ\beta\) operator described in [TS07]. This discontinuity capturing scheme involves the introduction of a dissipation term of the form

\[ \int_{\Omega} \nu_{SHOCK} \nabla \bm v \!:\! \nabla \bm q dV \]

The shock capturing viscosity is implemented following the first formulation described in {cite} tezduyar2007yzb. The characteristic velocity \(u_{cha}\) is taken to be the acoustic speed while the reference density \(\rho_{ref}\) is just the local density. Shock capturing viscosity is defined by the following

\[ \nu_{SHOCK} = \tau_{SHOCK} u_{cha}^2 \]

where,

\[ \tau_{SHOCK} = \frac{h_{SHOCK}}{2u_{cha}} \left( \frac{ \,|\, \nabla \rho \,|\, h_{SHOCK}}{\rho_{ref}} \right)^{\beta} \]

\(\beta\) is a tuning parameter set between 1 (smoother shocks) and 2 (sharper shocks. The parameter \(h_{SHOCK}\) is a length scale that is proportional to the element length in the direction of the density gradient unit vector. This density gradient unit vector is defined as \(\hat{\bm j} = \frac{\nabla \rho}{|\nabla \rho|}\). The original formulation of Tezduyar and Senga relies on the shape function gradient to define the element length scale, but this gradient is not available to qFunctions in libCEED. To avoid this problem, \(h_{SHOCK}\) is defined in the current implementation as

\[ h_{SHOCK} = 2 \left( C_{YZB} \,|\, \bm p \,|\, \right)^{-1} \]

where

\[ p_k = \hat{j}_i \frac{\partial \xi_i}{x_k} \]

The constant \(C_{YZB}\) is set to 0.1 for piecewise linear elements in the current implementation. Larger values approaching unity are expected with more robust stabilization and implicit timestepping.

Density Current

For this test problem (from [SWW+93]), we solve the full Navier-Stokes equations (15), for which a cold air bubble (of radius \(r_c\)) drops by convection in a neutrally stratified atmosphere. Its initial condition is defined in terms of the Exner pressure, \(\pi(\bm{x},t)\), and potential temperature, \(\theta(\bm{x},t)\), that relate to the state variables via

\[ \begin{aligned} \rho &= \frac{P_0}{( c_p - c_v)\theta(\bm{x},t)} \pi(\bm{x},t)^{\frac{c_v}{ c_p - c_v}} \, , \\ e &= c_v \theta(\bm{x},t) \pi(\bm{x},t) + \bm{u}\cdot \bm{u} /2 + g z \, , \end{aligned} \]

where \(P_0\) is the atmospheric pressure. For this problem, we have used no-slip and non-penetration boundary conditions for \(\bm{u}\), and no-flux for mass and energy densities.

Channel

A compressible channel flow. Analytical solution given in [Whi99]:

\[ u_1 = u_{\max} \left [ 1 - \left ( \frac{x_2}{H}\right)^2 \right] \quad \quad u_2 = u_3 = 0\]
\[T = T_w \left [ 1 + \frac{Pr \hat{E}c}{3} \left \{1 - \left(\frac{x_2}{H}\right)^4 \right \} \right]\]
\[p = p_0 - \frac{2\rho_0 u_{\max}^2 x_1}{Re_H H}\]

where \(H\) is the channel half-height, \(u_{\max}\) is the center velocity, \(T_w\) is the temperature at the wall, \(Pr=\frac{\mu}{c_p \kappa}\) is the Prandlt number, \(\hat E_c = \frac{u_{\max}^2}{c_p T_w}\) is the modified Eckert number, and \(Re_h = \frac{u_{\max}H}{\nu}\) is the Reynolds number.

Boundary conditions are periodic in the streamwise direction, and no-slip and non-penetration boundary conditions at the walls. The flow is driven by a body force determined analytically from the fluid properties and setup parameters \(H\) and \(u_{\max}\).

Flat Plate Boundary Layer

Laminar Boundary Layer - Blasius

Simulation of a laminar boundary layer flow, with the inflow being prescribed by a Blasius similarity solution. At the inflow, the velocity is prescribed by the Blasius soution profile, density is set constant, and temperature is allowed to float. Using weakT: true, density is allowed to float and temperature is set constant. At the outlet, a user-set pressure is used for pressure in the inviscid flux terms (all other inviscid flux terms use interior solution values). The wall is a no-slip, no-penetration, no-heat flux condition. The top of the domain is treated as an outflow and is tilted at a downward angle to ensure that flow is always exiting it.

Turbulent Boundary Layer

Simulating a turbulent boundary layer without modeling the turbulence requires resolving the turbulent flow structures. These structures may be introduced into the simulations either by allowing a laminar boundary layer naturally transition to turbulence, or imposing turbulent structures at the inflow. The latter approach has been taken here, specifically using a synthetic turbulence generation (STG) method.

Synthetic Turbulence Generation (STG) Boundary Condition

We use the STG method described in [SSST14]. Below follows a re-description of the formulation to match the present notation, and then a description of the implementation and usage.

Equation Formulation
\[ \bm{u}(\bm{x}, t) = \bm{\overline{u}}(\bm{x}) + \bm{C}(\bm{x}) \cdot \bm{v}' \]
\[ \begin{aligned} \bm{v}' &= 2 \sqrt{3/2} \sum^N_{n=1} \sqrt{q^n(\bm{x})} \bm{\sigma}^n \cos(\kappa^n \bm{d}^n \cdot \bm{\hat{x}}^n(\bm{x}, t) + \phi^n ) \\ \bm{\hat{x}}^n &= \left[(x - U_0 t)\max(2\kappa_{\min}/\kappa^n, 0.1) , y, z \right]^T \end{aligned} \]

Here, we define the number of wavemodes \(N\), set of random numbers \( \{\bm{\sigma}^n, \bm{d}^n, \phi^n\}_{n=1}^N\), the Cholesky decomposition of the Reynolds stress tensor \(\bm{C}\) (such that \(\bm{R} = \bm{CC}^T\) ), bulk velocity \(U_0\), wavemode amplitude \(q^n\), wavemode frequency \(\kappa^n\), and \(\kappa_{\min} = 0.5 \min_{\bm{x}} (\kappa_e)\).

\[ \kappa_e = \frac{2\pi}{\min(2d_w, 3.0 l_t)} \]

where \(l_t\) is the turbulence length scale, and \(d_w\) is the distance to the nearest wall.

The set of wavemode frequencies is defined by a geometric distribution:

\[ \kappa^n = \kappa_{\min} (1 + \alpha)^{n-1} \ , \quad \forall n=1, 2, ... , N \]

The wavemode amplitudes \(q^n\) are defined by a model energy spectrum \(E(\kappa)\):

\[ q^n = \frac{E(\kappa^n) \Delta \kappa^n}{\sum^N_{n=1} E(\kappa^n)\Delta \kappa^n} \ ,\quad \Delta \kappa^n = \kappa^n - \kappa^{n-1} \]
\[ E(\kappa) = \frac{(\kappa/\kappa_e)^4}{[1 + 2.4(\kappa/\kappa_e)^2]^{17/6}} f_\eta f_{\mathrm{cut}} \]
\[ f_\eta = \exp \left[-(12\kappa /\kappa_\eta)^2 \right], \quad f_\mathrm{cut} = \exp \left( - \left [ \frac{4\max(\kappa-0.9\kappa_\mathrm{cut}, 0)}{\kappa_\mathrm{cut}} \right]^3 \right) \]

\(\kappa_\eta\) represents turbulent dissipation frequency, and is given as \(2\pi (\nu^3/\varepsilon)^{-1/4}\) with \(\nu\) the kinematic viscosity and \(\varepsilon\) the turbulent dissipation. \(\kappa_\mathrm{cut}\) approximates the effective cutoff frequency of the mesh (viewing the mesh as a filter on solution over \(\Omega\)) and is given by:

\[ \kappa_\mathrm{cut} = \frac{2\pi}{ 2\min\{ [\max(h_y, h_z, 0.3h_{\max}) + 0.1 d_w], h_{\max} \} } \]

The enforcement of the boundary condition is identical to the blasius inflow; it weakly enforces velocity, with the option of weakly enforcing either density or temperature using the the -weakT flag.

Initialization Data Flow

Data flow for initializing function (which creates the context data struct) is given below:

flowchart LR subgraph STGInflow.dat y lt[l_t] eps Rij[R_ij] ubar end subgraph STGRand.dat rand[RN Set]; end subgraph User Input u0[U0]; end subgraph init[Create Context Function] ke[k_e] N; end lt --Calc-->ke --Calc-->kn y --Calc-->ke subgraph context[Context Data] yC[y] randC[RN Set] Cij[C_ij] u0 --Copy--> u0C[U0] kn[k^n]; ubarC[ubar] ltC[l_t] epsC[eps] end ubar --Copy--> ubarC; y --Copy--> yC; lt --Copy--> ltC; eps --Copy--> epsC; rand --Copy--> randC; rand --> N --Calc--> kn; Rij --Calc--> Cij[C_ij]

This is done once at runtime. The spatially-varying terms are then evaluated at each quadrature point on-the-fly, either by interpolation (for \(l_t\), \(\varepsilon\), \(C_{ij}\), and \(\overline{\bm u}\)) or by calculation (for \(q^n\)).

The STGInflow.dat file is a table of values at given distances from the wall. These values are then interpolated to a physical location (node or quadrature point). It has the following format:

[Total number of locations] 14
[d_w] [u_1] [u_2] [u_3] [R_11] [R_22] [R_33] [R_12] [R_13] [R_23] [sclr_1] [sclr_2] [l_t] [eps]

where each [  ] item is a number in scientific notation (ie. 3.1415E0), and sclr_1 and sclr_2 are reserved for turbulence modeling variables. They are not used in this example.

The STGRand.dat file is the table of the random number set, \(\{\bm{\sigma}^n, \bm{d}^n, \phi^n\}_{n=1}^N\). It has the format:

[Number of wavemodes] 7
[d_1] [d_2] [d_3] [phi] [sigma_1] [sigma_2] [sigma_3]

The following table is presented to help clarify the dimensionality of the numerous terms in the STG formulation.

Math

Label

\(f(\bm{x})\)?

\(f(n)\)?

\( \{\bm{\sigma}^n, \bm{d}^n, \phi^n\}_{n=1}^N\)

RN Set

No

Yes

\(\bm{\overline{u}}\)

ubar

Yes

No

\(U_0\)

U0

No

No

\(l_t\)

l_t

Yes

No

\(\varepsilon\)

eps

Yes

No

\(\bm{R}\)

R_ij

Yes

No

\(\bm{C}\)

C_ij

Yes

No

\(q^n\)

q^n

Yes

Yes

\(\{\kappa^n\}_{n=1}^N\)

k^n

No

Yes

\(h_i\)

h_i

Yes

No

\(d_w\)

d_w

Yes

No

Meshing

The flat plate boundary layer example has custom meshing features to better resolve the flow. One of those is tilting the top of the domain, allowing for it to be a outflow boundary condition. The angle of this tilt is controled by -platemesh_top_angle

The primary meshing feature is the ability to grade the mesh, providing better resolution near the wall. There are two methods to do this; algorithmically, or specifying the node locations via a file. Algorithmically, a base node distribution is defined at the inlet (assumed to be \(\min(x)\)) and then linearly stretched/squeezed to match the slanted top boundary condition. Nodes are placed such that -platemesh_Ndelta elements are within -platemesh_refine_height of the wall. They are placed such that the element height matches a geometric growth ratio defined by -platemesh_growth. The remaining elements are then distributed from -platemesh_refine_height to the top of the domain linearly in logarithmic space.

Alternatively, a file may be specified containing the locations of each node. The file should be newline delimited, with the first line specifying the number of points and the rest being the locations of the nodes. The node locations used exactly at the inlet (assumed to be \(\min(x)\)) and linearly stretched/squeezed to match the slanted top boundary condition. The file is specified via -platemesh_y_node_locs_path. If this flag is given an empty string, then the algorithmic approach will be performed.