xref: /libCEED/examples/fluids/README.md (revision 4de8550ab4506b5c5bd81d0c3febcb71b9614dca)

## libCEED: Navier-Stokes Example

This page provides a description of the Navier-Stokes example for the libCEED library, based on PETSc.
PETSc v3.17 or a development version of PETSc at commit 0e95d842 or later is required.

The Navier-Stokes problem solves the compressible Navier-Stokes equations in three dimensions using an explicit time integration.
The state variables are mass density, momentum density, and energy density.

The main Navier-Stokes solver for libCEED is defined in [`navierstokes.c`](navierstokes.c) with different problem definitions according to the application of interest.

Build by using:

`make`

and run with:

```
./navierstokes -ceed [ceed] -problem [problem type] -degree [degree]
```

## Runtime options

% inclusion-fluids-marker

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

:::{list-table} Common Runtime Options
:header-rows: 1

* - 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
  - `0`

* - `-ts_monitor_solution`
  - PETSc output format, such as `cgns:output-%d.cgns` (requires PETSc `--download-cgns`)
  -

* - `-ts_monitor_solution_interval`
  - Number of time steps between visualization output frames.
  - `1`

* - `-viewer_cgns_batch_size`
  - Number of frames written per CGNS file if the CGNS file name includes a format specifier (`%d`).
  - `20`

* - `-checkpoint_interval`
  - Number of steps between writing binary checkpoints. `0` has no output, `-1` outputs final state only
  - `10`

* - `-checkpoint_vtk`
  - Checkpoints include VTK (`*.vtu`) files for visualization. Consider `-ts_monitor_solution`instead.
  - `false`

* - `-viz_refine`
  - Use regular refinement for VTK visualization
  - `0`

* - `-output_dir`
  - Output directory for binary checkpoints and VTK files (if enabled).
  - `.`

* - `-output_add_stepnum2bin`
  - Whether to add step numbers to output binary files
  - `false`

* - `-continue`
  - Continue from previous solution (input is step number of previous solution)
  - `0`

* - `-continue_filename`
  - Path to solution binary file from which to continue from
  - `[output_dir]/ns-solution.bin`

* - `-continue_time_filename`
  - Path to time stamp binary file (only for legacy checkpoints)
  - `[output_dir]/ns-time.bin`

* - `-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
  -

* - `-bc_freestream`
  - Use freestream 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`, `bc_freestream`  and/or `-bc_slip_x`, `-bc_slip_y`, and `-bc_slip_z` are:

:::{list-table} 2D Face ID Labels
:header-rows: 1
* - PETSc Face Name
  - Cartesian direction
  - Face ID

* - faceMarkerBottom
  - -z
  - 1

* - faceMarkerRight
  - +x
  - 2

* - faceMarkerTop
  - +z
  - 3

* - faceMarkerLeft
  - -x
  - 4
:::

:::{list-table} 2D Face ID Labels
:header-rows: 1
* - 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:

:::{list-table} Advection2D Runtime Options
:header-rows: 1

* - 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:

:::{list-table} Advection3D Runtime Options
:header-rows: 1

* - 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:

:::{list-table} Isentropic Vortex Runtime Options
:header-rows: 1

* - 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:

:::{list-table} Shock Tube Runtime Options
:header-rows: 1

* - 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:

:::{list-table} Newtonian Ideal Gas problems Runtime Options
:header-rows: 1

* - 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, 60, 128 for degree = 1, 2, 3`
  -

* - `-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

* - `-state_var`
  - State variables to solve solution with. `conservative` ($\rho, \rho \bm{u}, \rho e$) or `primitive` ($P, \bm{u}, T$)
  - `conservative`
  - string
:::

#### Newtonian Wave

The newtonian wave problem has the following command-line options in addition to the Newtonian Ideal Gas options:

:::{list-table} Newtonian Wave Runtime Options
:header-rows: 1

* - Option
  - Description
  - Default value
  - Unit

* - `-freestream_riemann`
  - Riemann solver for boundaries (HLL or HLLC)
  - `hllc`
  -

* - `-freestream_velocity`
  - Freestream velocity vector
  - `0,0,0`
  - `m/s`

* - `-freestream_temperature`
  - Freestream temperature
  - `288`
  - `K`

* - `-freestream_pressure`
  - Freestream pressure
  - `1.01e5`
  - `Pa`

* - `-epicenter`
  - Coordinates of center of perturbation
  - `0,0,0`
  - `m`

* - `-amplitude`
  - Amplitude of the perturbation
  - `0.1`
  -

* - `-width`
  - Width parameter of the perturbation
  - `0.002`
  - `m`

:::

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

```
./navierstokes -options_file newtonianwave.yaml
```

```{literalinclude} ../../../../../examples/fluids/newtonianwave.yaml
:language: yaml
```

#### Vortex Shedding - Flow past Cylinder

The vortex shedding, flow past cylinder problem has the following command-line options in addition to the Newtonian Ideal Gas options:

:::{list-table} Vortex Shedding Runtime Options
:header-rows: 1

* - Option
  - Description
  - Default value
  - Unit

* - `-freestream_velocity`
  - Freestream velocity vector
  - `0,0,0`
  - `m/s`

* - `-freestream_temperature`
  - Freestream temperature
  - `288`
  - `K`

* - `-freestream_pressure`
  - Freestream pressure
  - `1.01e5`
  - `Pa`

:::

The initial condition is taken from `-reference_temperature` and `-reference_pressure`.
To run this problem, first generate a mesh:

```console
$ make -C examples/fluids/meshes
```

Then run by building the executable and running:

```console
$ make build/fluids-navierstokes
$ mpiexec -n 6 build/fluids-navierstokes -options_file vortexshedding.yaml
```

The vortex shedding period is roughly 6 and this problem runs until time 100 (2000 time steps).

```{literalinclude} ../../../../../examples/fluids/vortexshedding.yaml
:language: yaml
```

#### Density current

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

:::{list-table} Density Current Runtime Options
:header-rows: 1

* - 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 has the following command-line options in addition to the Newtonian Ideal Gas options:

:::{list-table} Channel Runtime Options
:header-rows: 1

* - 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
```
```{literalinclude} ../../../../../examples/fluids/channel.yaml
:language: yaml
```

#### Blasius boundary layer

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

:::{list-table} Blasius Runtime Options
:header-rows: 1

* - Option
  - Description
  - Default value
  - Unit

* - `-velocity_infinity`
  - Freestream velocity
  - `40`
  - `m/s`

* - `-temperature_infinity`
  - Freestream temperature
  - `288`
  - `K`

* - `-temperature_wall`
  - Wall temperature
  - `288`
  - `K`

* - `-delta0`
  - Boundary layer height at the inflow
  - `4.2e-3`
  - `m`

* - `-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.
  - `""`
  -

* - `-n_chebyshev`
  - Number of Chebyshev terms
  - `20`
  -

* - `-chebyshev_`
  - Prefix for Chebyshev snes solve
  -
  -

:::

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

```
./navierstokes -options_file blasius.yaml
```

```{literalinclude} ../../../../../examples/fluids/blasius.yaml
:language: yaml
```

#### STG Inflow for Flat Plate

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

:::{list-table} Blasius Runtime Options
:header-rows: 1

* - 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`
  -

* - `-stg_fluctuating_IC`
  - "Extrude" the fluctuations through the domain as an initial 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.