Unstructured Solver Inputs

Input File Syntax

LAVA utilizes a plain-text input file to specify solver options. The input file is organized in several sections.

Info

All the inputs are interpreted in standard SI units

Most inputs have default values built into the code and can be skipped from the input file if default value is acceptable.

The input file is case-insensitive, with the exception of string entries such as file names. In addition, most white spaces are disregarded with the exception of carriage return.

Input entries are separated by either a semicolon ";" or carriage return. To illustrate, let's look at two examples:

ExampleInput1 = 1.0
ExampleInput2 = 2.0

is equivalent to

example input 2 = 2.0; example input 1 = 1.0

Note that the order of entries are swapped in the second example.

The file is organized in nested sections that encapsulate logical groups of inputs together. The sections are denoted with a section name followed by an opening curly brace "{" and ending with a closing curly brace "}". Here is an example:

Section1 {
    input1 = 1
    SubsectionA {foo = 2}
    SubsectionB {bar = 3}
}

This example shows a section and has one input entry as well as two nested sections both of which have their own inputs entries.

Like ordinary input entries, the order of sections in the input file is not important as long as the nesting hierarchy is correct.

The last piece of relevant syntax is the numbered sections that are utilized to construct numbered lists, suffixed by an underscore "_" followed by item number. One example usage is listing of boundary conditions for the unstructured solver:

BC_1 {type=wall}
BC_2 {type=inlet}
BC_3 {type=farfield}

Vector entries (or lists) are denoted using brackets:

velocity = [2,5,3]

Lists follow the same format as vectors:

velocity = [2,3,5-11,14]

Note that you can specify ranges (inclusive) with dashed as above.

Comment lines can be added to the input file using hash signs # either in a separate line or in-line with the code:

# This is a comment
entry = 3 # this is also a comment

Fluid Section

The fluid, either a single species, or a mixture of species is specified in the fluid section:

fluid {
    species = air_perfect
}

This example uses perfect air as an equivalent single-species gas. The code is distributed with a database of common gas species. The database is stored in individual "species files" that contain gas properties such as molecular weight, specific heat ratio, viscosity (or Sutherland coefficients), Prandtl number etc...

Available fluid species are: air_perfect, h (hydrogen), o (oxygen), h2, o2, oh,h2o (water vapor),co and co2. Users can create their own species definitions.

For multi-species problems, you can specify the mixture composition instead of single species such as:

fluid {
    mixture = [h,o,h2,o2,h2o,oh]
}

Info

LAVA currently does not support reacting flows.

Reference Conditions Section

Reference conditions list the reference flow quantities (typically the freestream conditions in external flow problems) as well as the reference area, length and angles (pitch and yaw) in airplane coordinates. These values are used in initial and boundary conditions (unless different conditions are explicitly specified). They are also used to normalize force and moment coefficients.

The reference conditions need to properly specify the thermodynamic state and velocity field. Typically two of the three thermodynamic quantities (density, pressure and temperature) need to be specified. Specifying all three of them will result in an error message as it is an over-specification. In addition, velocity vector or speed and angles are typically specified. Here are a few example ways to specify the reference conditions:

reference conditions {
    pressure    = 1.e5
    temperature = 300
    velocity    = [100.,0.,0.]
    area        = 1.0
    length      = 1.0
    # If using multi-species, list the mass fractions of each species
    mass fractions = [1,0,0]
}

Another way to specify the thermodynamic conditions is via prescribe an altitude in the standard atmosphere model. Using this requires the flag use std atmosphere = true.

reference conditions {
    use standard atmosphere = true # false by default
    std altitude = 1000
    velocity     = [100.,0.,0.]
    area         = 1.0
    length       = 1.0
    # If using multi-species, list the mass fractions of each species
    mass fractions = [1,0,0]
}

Info

The thermodynamic condition description for multi-phase is slightly different than for single and multi-species. Please refer to the Cartesian Multiphase Solver page for details.

By default, reference area and length values are unity (1.0) and can be skipped if normalized forces are not needed. Instead of the velocity, we can utilize Mach number and flow angles:

reference conditions {
    pressure    = 1.e5
    temperature = 300
    Mach        = 0.5
    alpha       = 4
    beta        = 0
}

Note that alpha and beta are pitch and yaw angles in degrees. This assumes that the model is in the airplane coordinate frame (+x: from nose to tail, +y: along the right wing, +z: from airplane towards the ground). These angles do not transform the model, but rather set the freestream direction so that the model is at the desired angles relative to the freestream. Positive alpha corresponds to pitch up, while positive beta corresponds to yaw towards right. The following coordinate transform is used in the code to transform the reference velocity vector to a Cartesian frame.

\[ \begin{align*} v_x &=v_0 \cdot \cos(\alpha) \cdot \cos(\beta)\\ v_y &= -v_0 \cdot \sin(\beta)\\ v_z &=v_0 \cdot \sin(\alpha) \cdot \cos(\beta) \end{align*} \]

For launch-vehicle simulations, the flow angle can be specified via alpha, angle of attack, and phi, the roll angle. This assumes that the rocket model is aligned in the x axis, this +x being from windward to leeward. phi=0 corresponds to incoming flow aligned to the x-z plane. alpha then indicates the rotation of the incoming flow on the -y axis and phi the rotation of the incoming flow on the -x axis. The following coordinate transform is used in the code to transform the reference velocity vector to a Cartesian frame.

\[ \begin{align*} v_x &=v_0 \cdot \cos(\alpha)\\ v_y &=v_0 \cdot \sin(\alpha)\cdot\sin(\phi) \\ v_z &=v_0 \cdot \sin(\alpha) \cdot \cos(\phi) \end{align*} \]

Another way to specify the reference conditions is by using Reynolds number, Mach number and temperature:

reference conditions {
    Mach        = 1.4
    Re          = 5.e6
    alpha       = 4
    beta        = 2
    length      = 1.8
}

Note that the length is the reference length used in the Reynolds number definition.

Yet another way to specify the reference conditions is by using dynamic pressure, Mach number, and standard atmosphere. Using this requires the flag use std atmosphere = true.

reference conditions {
    use std atmosphere = true
    dynamic pressure   = 70000
    Mach        = 1.4
    alpha       = 4
    beta        = 2
}

When the length scale and area are specified, the force and moment coefficients (if loads output is requested), will be utilized these refence values for normalization. For external airplane aerodynamics, for example, length would correspond to the mean chord and area would correspond to the main wing planform area.

Warning

Correct rolling moment coefficient should use airplane span length for normalization. However, only one length scale input is currently supported and the rolling moment coefficient that LAVA outputs will also be normalized by the chord length. It is up to the user to make the rolling moment coefficient correction.

reference conditions {
    ...
    area        = 280.
    length      = 5.
    dx          = 1.e-3
}

The dx input represents a characteristic grid spacing value, often indicating the cell size (in meters) in the finest region of the grid. It is not required (default value is 1), but it can make life easier when setting up some of the other input sections (e.g. trip lines, relative grid spacing w.r.t. dx). The role of dx will be described in relevant sections that make use of this parameter.

For multi-species runs, species mass fractions should be specified:

reference conditions {
    ...
    mass fractions = [1,0,0]
}

Time Section

The time section specifies whether the run is steady or unsteady, and if unsteady, sets up the extents of step size for the time integration.

In the simplest case of steady runs, the section can be skipped all together or the section can be used just to explicitly state a steady setup:

time {
    steady = true
}

For unsteady runs, several other entries become revelant. The following example setups up a case that will run for 1000 time steps with a constant dt.

time {
  steady = false
  number of steps = 1000
  dt = 1.e-5
}

Alternatively, you can specify stop condition in terms of total time integrated. You can also choose to start the time integration from a non-zero inital time:

time {
  steady = false
  start time  = 0.1
  finish time = 1.
  dt = 1.e-5
}

Yet another convenient and common way to specify time integration bounds is by utilizing CTUs (Convective Time Unit). A CTU is simply the non-dimensional time normalized by reference flow velocity magnitue and the reference length as specified in the Reference Conditions Section.

time {
  steady = false
  start ctu  = 1.
  finish ctu = 10.
  dt = 1.e-5
}

Time step size can be specified in terms of CFL instead of a fixed dt as used above. The specifed cfl number is interpreted as the maximum allowed CFL in the whole domain and the corresponding dt value is adjusted automatically.

time {
  steady = false
  finish ctu  = 90
  initial cfl = 0.05
  cfl = 2.
}

Here initial cfl value is used in the start-up and cfl gradually climbs until the specified cfl number is met.

Unstructured Specific Inputs

Unstructured specific inputs are grouped within the unst {} section. Every following section here should be included within it.

Grid File and Transformations

The grid is expected to be in CGNS format, written as arbitrary polyhedrals (even if the grid is made up of only standard element types). Only a single volume zone is supported while each following surface zone is interpreted as a separate boundary condition (BC).

The grid file name is specified as:

grid {
    file = grid.cgns
}

Note that the specified grid file (or a symbolic link to it) needs to be present in the same directory as the input file. Otherwise, you can enter the full file path.

Multiple grid files can be specified. The solver will automatically merge the multiple grids and associate each with a domain id. This is typically used in sliding mesh simulations with multiple moving/stationary domains.

grid {
    file = [grid1.cgns, grid2.cgns, grid3.cgns]
}

Some basic grid transformation operations are possible within the grid file. See the examples below:

grid {
    file = grid.cgns
    # Scale the grid by a factor of 2
    transform_1 {function=scale; factor=2}
    # Translate the grid 1 unit length in x direction
    transform_2 {function=translate; distance=[1,0,0]}
    # Rotate the grid 90 degrees around y axis with an origin of (0,0,0)
    transform_3 {function=rotate; center=[0,0,0]; axis=[0,1,0]; angle=90}
    # Convert 1-cell thick (essentially 2D) grid into axisymmetric (wedge shape)
    # Depth is assumed to be in the depth direction, while y-coordinate is the height.
    transform_4 {function=axisymmetric; y_center=0, y_tolerance=1.e-6, angle=5}
}

Info

Everything else in the input file and the case setup need to be with respect to the already transformed grid.

The transformations are applied in the given order. You can list as many transformations as needed. If no transformations are needed, the transform inputs can be skipped.

The unstructured solver operates on "tiles" of cells for computational efficiency. Size of these tiles, i.e. approximately how many cells there should be in each tile, can be specified as:

grid {
    file = grid.cgns
    tile size = 1024
}

The tile size affects the computational performance of the code, and its optimum value depends on the type of grid and the architecture of the system being used. There are built-in default "best practices" for several NAS systems so you may choose to skip specifying a tile size in favor of the defaults. In general, especially if you are running large scale cases, it is best to experimentally determine the best choice for tile size.

Additional partitioning controls are available to customize the METIS partitioning algorithm used for tile creation:

grid {
    file = grid.cgns
    tile size = 1024
    # Enable cell-based weighting during partitioning
    # When true, uses number of faces per cell as weight for load balancing
    # When false, treats all cells with equal weight (default: false)
    do partition weights = false

    # Choose partitioning method for tile creation
    # true: Use graph-based partitioning (cell-to-cell connectivity)
    # false: Use dual mesh partitioning (face-based connectivity) (default: false)
    do graph tile partition = false
}

These partitioning options are optional and control how METIS creates computational tiles. The do partition weights option defaults to false, meaning all cells are treated with equal weight during partitioning. When enabled, it uses the number of faces per cell as a weight, which can improve load balancing by accounting for computational cost differences between cell types (e.g., tetrahedra vs hexahedra). The do graph tile partition option also defaults to false, using dual mesh partitioning where faces are treated as nodes. When set to true, it switches to graph-based partitioning using direct cell-to-cell adjacency, which typically produces higher quality partitions but may have different performance characteristics depending on mesh topology.

Note

If the unstructured mesh in question was a Voronoi tessellation of the computational domain created by the LAVA Voronoi Mesher (see the Voronoi Mesh Generator section for details), then the coordinates of the cell generating points, i.e. seeds, are provided in the mesh file. In this case, the user has the choice to read the seed coordinates information from the *.cgns file, the choice to use these coordinates instead of the cell centroids and the choice to use the location where the line segment between two neighboring seeds intersect their respective face instead of the face centroid as the flux exchange location. Additionally, the user can choose to include a correction to green-gauss gradient operator exclusive to Voronoi-type meshes. This correction was explained by Volker Springel in his 2009 paper. These behaviors are controlled by the following parameters, all of which are disabled by default.

grid {
    read seed coords = false #default
    use seed coords = false #default
    update face centroids = false #default
    do voro grad correction = false #default
}

By default, grid diagnostics is active (see the Diagnostics section for details).

grid {
    do diagnostics = true # default
}

Diagnostics

At runtime, solver will analyze the grid to determine problem cells. In most cases, these are cell that are deemed degenerate or problematic for solver stability. The cells marked as problem are filtered out, meaning they are not solved for. Instead, they are updated at each step by averaging their neighbor values. A summary of how many cells were filtered is provided at run time.

The diagnostics metrics and default threshold values are specified as listed below. Note that the defaults are designed to catch the most severe problems and have been calibrated to ensure solver robustness. In most cases, keeping the default values are recommended. If you are using defaults, this section or any default entry is not required to be present in the input file.

diagnostics {
  # Maximum allowed angle (in degrees) between face normal and
  # the vector connecting the two adjacent cell centers
  max nonorthogonality angle = 90.
  # Maximum order of magnitue change of volume between adjacent cells
  max volume change order  = 4.
  # Checks accuracy of cell gradient stencil/weights vs a linear function
  # The error is normalized and specified as percentage.
  max gradient error       = 1.e3
  # Check convexity of cells
  do convexity check       = true
  # Check if a cell is watertight
  max green gauss error    = 1e-8 # 1e-3 for FLOAT executable
  # Check if all faces have a valid neighbor cell
  # if the face is not a BC and a neighbor face was not found, abort the run
  do face neighbor check = true;
  # Check how skewed the angle from cell-centroid-to-neighbor
  # versus the cell-centroid-to-face-centroid vector is. Input in degrees.
  do neighbor face angle check = false
  min face centroid angle = 45
  # Check how close the cell centroid is to its faces. Diagnostic is normalized
  # as (min face dist)/(max face dist) for a pair of cells on either side
  # of a face.
  do face dist growth check = false
  min face dist growth = 4
  # Check isotropy of the cell. Closer to 1 means perfectly isotropic.
  do iso factor check = false
  min iso factor = 0.25
  # If the following is set as nonzero, then cell diagnotics are tracked into
  # histograms with the given number of bins  and output as CSV files.
  # There is one CSV file per test, stored in
  # the folder "diagnostics" under the file names
  #
  # face_planarity.csv, face_orthogonality.csv, cell_convexity.csv,
  # face_alignment.csv, cell_size_filtering.csv, volumetric_jump.csv,
  # green_gauss.csv, gradient_err.csv,  face_neighbor_angle.csv,
  # face_equidistant.csv, and cell_isotropy.csv.
  #
  # The green_gauss.csv is given in log scale. The cell_size_filtering.csv
  # always uses 10% filter size bin widths, matching the terminal output.
  histogram bins = 0
  # If the following is activated, user manually mark additional problem cells
  # by listing their global ids (ids as in the provided grid file) in an ASCII file
  #problem cell file = problemcells.txt
  #
  # When enabled, write the geometry and flow solution of any problem cells to
  # volume/Unst.problemCells.<step>.cgns at every startup (fresh start and restart).
  # On an abort triggered by robust mode, the file is named
  # volume/ABORT_Unst.problemCells.<step>.cgns instead.
  # Output is skipped automatically if no problem cells exist.
  # This option is TRUE by default
  dump problem cells = false
}

Tip

At every startup, if any problem cells are found, the solver automatically writes their geometry and flow solution to volume/Unst.problemCells.<step>.cgns. This file can be opened in ParaView to visually inspect where the problematic cells are located in the mesh.

Warning

If you encounter a segfault during problem cell output, please try tunring it off and running again.

Another use of the diagnostics for runs using explicit time integration is to filter-out cells that are too small. The filtered cells do not contribute to CFL computation and do not cause time step size restrictions. You can leave this unspecified if not using explicit time integration.

diagnostics {
    # Specify minimum cell size intended when grid was created
    target size filter = 4e-3
    # Solver will use a histogram approach and determine the minimum
    # allowed cell size based on maximum percentage of cells to be filtered due to size
    max percent cells filtered = 0.05 # default
    # You can also manually specify a hard cut-off for the
    # minimum allowed cell size. By default, it is half of the target size filter
    min cell size filter = 2e-3
}

Robust Mode

Note

Some features of the robust mode is currently only supported for explicit time integration

In rare cases where solver is unstable, you can turn on the robust mode. In this mode, at every time integration step (or Runge-Kutta integrator stage), the solver will test the solution update against the Mach, density, pressure and temperature limits as specified in the navier stokes section. If any of these criteria are violated, the solver will attempt to recover by applying a smoothing operator to the problematic cells and repeating the step by using dissipative numerics for those cells for the repeated step only.

expert {
    robust mode = true # false by default
    max bomb recovery attempts = 5
    # Permanently add the crashing cell to the problem cell list
    mark bomb cells as problem = false
    # When enabled, write the geometry and flow solution of the offending tiles to
    # volume/bombTiles.cgns on abort
    # This option is TRUE by default
    dump bomb tiles = false
}

If the stability issue persists after max bomb recovery attempts consecutive attempts, the solver will abort the solution and write out the problematic tile solutions to volume/bombTiles.cgns file.

If mark bomb cells as problem = true, the cells experience the stability issue are permanently added to the problem list and are no longer solved for.

Tip

On abort, the solver writes volume/ABORT_Unst.problemCells.<step>.cgns containing the flow solution at the time of abort. If mark bomb cells as problem = true is set in the input file, the cell that caused the crash will also be a part of the problem cell list. Open this file in ParaView and use the Process Unst Mask filter from the LAVAParaviewPlugin (please refer to the plugin documentation here) to extract more details about the cell.

Warning

In implicit simulations, robust mode only considers the maximum mach specification in the navier stokes section. mark bomb cells as problem and max bomb recovery attempts inputs are not used.

Periodic

The domain can be setup as periodic in x,y,z directions. However the grid must be built with special considerations in one of two ways. One option is to apply extra grid cell padding (at least two cells thick on both periodic ends) to act as periodic fringe exchange locations. Ideally, these padding cells, when shifted in the periodic direction by the domain width, should match interior cells perfectly on the other end. There are tools to mirror and extrude the cells to fit this criteria. Alternatively, the connectivity between cells at the periodic boundary must be correctly accounted for during the mesh generation step. This can be achieved in the LAVA Voronoi mesher, for example, by creating ghost periodic cells with the correct parent id where necessary. The end product is a seamless periodic domain generation.

If the first option is chosen and the grid is prepared, inputs for periodic setup are:

periodic {
    # Activate periodicity in y-direction
    do y = true
    # Interior domain bounds in periodic direction
    ymin = -0.0325
    ymax = 0.0325
    # In case periodic fringes do not exactly overlap with donor cells,
    # perform interpolation
    do interpolation = false # default
    # Test the periodic setup by performing an interpolation
    # for coordinates and checking linear exactness
    # false by default
    do test projection = true
}

The example above is only for y-direction. Periodicity can be activated in other directions at the same time. Here is an example for a triple periodic setup:

periodic {
    # Activate periodicity in all 3 directions
    do x = true
    do y = true
    do z = true
    # Interior domain bounds in periodic direction
    xmin = -0.5
    xmax = 0.5
    ymin = -0.5
    ymax = 0.5
    zmin = -0.5
    zmax = 0.5
}

Alternatively, if the periodicity was included in the mesh generation process in a way that correctly accounts for the periodicity and connectivity of the cells, then the user should inform the domain boundaries in each axial direction and turn on the face translation corrections via the following parameters:

periodic {
    do lava voro periodic = true

    # Interior domain bounds in periodic direction
    xmin = -0.5
    xmax = 0.5
    ymin = -0.5
    ymax = 0.5
    zmin = -0.5
    zmax = 0.5
}

Restarting

To resume runs, all you need to do is point to the restart file location as:

restart {
    file = ./restart/Unst_restartLatest.hdf5
}

The restart files (also referred to as checkpoint files) are saved periodically in restart directory created by the code under the run directory. The interval at which the restart files are created can be specified by users. The details are provided in restart output options section.

If the specified restart file is not found, the solver will print a message and proceed with a fresh start instead of aborting. This allows the same input file to be used for both initial runs and restarts — if the restart file does not yet exist, the simulation simply starts from scratch. If restart is not desired at all, you can skip or comment out this section.

Additional restart controls are:

restart {
    file = ./restart/Unst_restartLatest.hdf5
    # Set physical time to 0 on restart
    reset time = true
    # Set the time step number to 0 on restart
    reset nt   = true
}

Explicit Stepping

Explicit time stepping is available for unsteady simulations, typically used for WMLES or LES.

explicit {
    # Specify the type of Runge-Kutta time integration scheme
    # Options: rk1, rk4, ssprk22, ssprk32, ssprk33, ssprk53
    method = ssprk33 # default
}

Implicit Stepping

Implicit stepping options for either steady simulations, or pseudo-time stepping unsteady simulations are specified in the implicit section.

implicit {
    # method options are either 1st order accurate backward differencing formula (BDF1) or
    # the second order version BDF2. The default is BDF2
    # This option is not relevant for steady runs and can be skipped.
    method = BDF2
    # Number of implicit iterations.
    # For steady runs, this is the maximum number of iterations
    # For unsteady runs, this is the maximum number of sub-iterations for each time step
    number of steps = 10000
    # Specify CFL value
    CFL = 10
}

CFL entry is interpreted as local CFL, meaning that every cell has a different pseudo-time step size that is computed to match the CFL number in every cell.

In most cases, we favor using adaptive local CFL time stepping. This allows the code to try larger step sizes while adjusting as needed. Adaptive stepping generally improves robustness and results in faster convergence.

    implicit {
        ...
        # Allow CFL to adaptively vary between min and max
        adaptive {
            min CFL = 0.5
            max CFL = 100
            # In cases where turbulence model is used, CFL for
            # turbulence model equations can be limited differently
            max CFL turbulence = 50
        }
    }
}

Following are some options for the evaluation of Jacobian and the non-linear convergence

implicit {
    ...
    # Compute the jacobian analytically (approximate) or numerically (via finite different)
    jacobian method = analytical # analytical (default) or numerical
    # Interval of implicit steps at which the jacobian is updated
    jacobian interval = 5 # default
    # For numerical Jacobians, we can use the Q vector to automatically
    # determine the choice of epsilon. This input controls that option.
    do auto epsilon = true #default
    # If auto epsilon set false, then instead read an input epsilon value.
    jacobian epsilon  = 1.e-6 # default
}

By default, Jacobian method is analytical where analytical derivative of van Leer convective flux is used. This is approximate and evaluated without accounting for the second order MUSCL reconstruction. This is the most robust option although it can limit the depth of convergence in complicated cases. Numerical Jacobians are evaluated by finite differencing the RHS operator. While this method can be more accurate, it is more expensive and may be less robust due to finite differencing precision. The epsilon value is what is used to perturb primitive variables (normalized by reference quantities) when computing finite difference (numerical) Jacobians. It is computed automatically by default, but may be set to a fixed value if desired.

Finally, you may choose to activate a line search algorithm for the non-linear implicit iterations to find asymptotic convergence path.

implicit {
    ...
    # Whether to use line search
    do line search = false # default
    # Maximum relaxation factor for line search updates
    max relax = 1. # default
}

By default, the line search is turned-off as it is currently not reliable in general but it can improve the convergence rate if successful. max relax entry adjusts what portion of the line search update to advance the solution with.

Convergence Criteria

Convergence tolerances for implicit iterations can be specified as:

convergence {
    test continuity only = true #default
    relative tolerance = 1.e-6
    absolute tolerance = 1.e-12
}

The values are L2 error norms either relative to the first iteration (normalized), or in absolute (non-normalized) sense. The implicit iterations are deemed converged when all L2 error norms for all the equations meet either the relative or absolute criteria. If the test continuity only flag is set to true (default value), L2 error norm of only the continuity equation is considered as convergence criteria.

If these inputs are not provided, the code will simply perform the full number of steps specified in the implicit section regardless of convergence level.

Convective Terms

Convective flux options are specified in the convective terms section. Three different convective schemes can be setup; MUSCL, CENTRALQ and KEEPBLEND, latter two of which are low dissipation schemes to be used in scale-resolving simulations (i.e. LES and WMLES).

MUSCL Scheme

The traditional \(2^\mathrm{nd}\) order MUSCL scheme with limited gradient reconstruction can be setup with:

convective terms {
    # Convective flux scheme (options = MUSCL, CentralQ, KEEPBLEND) [default = MUSCL]
    scheme = MUSCL
    # MUSCL reconstruction order [options = 1,2 | default = 2]
    order = 2
    # Limiter [options = none, minmod, minmod soft | default = none]
    limiter = minmod
    # Only active for minmod soft, the overshoot parameter will allow normalized
    # overshoots of the specified amount. In this example, it is set to 5%
    limiter overshoot = 0.05
    # Flux function [options = roep,ausmpw+,hrslau2,rusanov | default = ausmpw+]
    flux function = ausmpw+
}

Choice of flux function can have a large impact in simulation accuracy and stability. Our current best practice is to use the ausmpw+ flux at all speeds. At high Mach numbers, we recommend adding extra dissipation (see High Mach Simulations section).

roep option stands for preconditioned Roe flux function, which has additional options:

convective terms {
    ...
    flux function = roep
    # Entropy fix method
    # 0: No entropy fix
    # 1: Entropy fix version 1
    # 2: Entropy fix version 2 (default)
    ientropy = 2   # default
    # Activate low-speed preconditioning
    # 1: No preconditioning (default)
    # 2: Local inviscid velocity
    # 3: Constant beta
    iprecon  = 0   # default
    # For iprecon = 2, beta is a measure of local Mach number (but w.r.t to the reference speed of sound)
    # Local speed of sound is scaled according to this Mach number
    beta     = 1.  # default
    # Epsilon threshold for eigen values in entropy fix
    eigeps   = 0.5 # default
}

limiter is only required for transonic or supersonic cases to preserve monotonicity around discontinuities.

Scale-Resolving Simulations

Scale-resolving simulations require reduced dissipation numerics to preserve smaller turbulent fluctuations that are at a comparable scale to individual grid cells. LAVA unstructured offers two different low dissipation schemes: CENTRALQ and KEEPBLEND. Both of them utilize blending sensors to selectively reduce dissipation in smooth flow regions (free of shocks) where turbulent fluctuations are present.

Upwind Blending Sensors

These sensors dictate where and to what degree the scheme deploys upwind dissipation. There are a variety of sensor options:

convective terms {
    # Sensor(s) that determine the amount of central Q blending or flux blending
    # Options: ducros, strain, limiter, ducrosmod, strainmod, grid quality
    # More than one sensor can be activated through the list below
    upwind blending sensors = [ducrosmod, limiter]
    # For ducrosmod and strainmod sensors, this parameter adjust the localization of the sensor
    # as a multiplier for the enstrophy term
    # Higher values represent more localization (i.e. less dissipation)
    # A value of 1 recovers original formulation, e.g. ducrosmod becomes the standard ducros sensor
    upwind sensor enstrophy factor = 100 # [default]
    # Limits of blending
    # 0: no upwind (full central), 1: full upwind
    upwind blending min          = 0. # default: 0
    upwind blending max          = 1. # default: 1
    # Same limits but only applied at wall adjacent faces
    wall upwind blending min     = 0. # default: 0
    wall upwind blending max     = 0. # default: 1
}

Our current best practice is to use the ducrosmod sensor with upwind sensor enstrophy factor = 100 for low Mach flows, and reduced factor values as needed for high Mach flows. At high supersonic/hypersonic regimes, upwind sensor enstrophy factor = 1 recovers the original Ducros sensor formulation and provides robustness.

CENTRALQ Scheme

This scheme blends left and right reconstructed primitive flow variables based on the upwind sensor. The details can be found in Sozer et al.¹. Example usage is shown below:

convective terms {
    scheme  = CentralQ
    order   = 2
    limiter = minmod
    # Pick which Q components to centralize (options: p,u,v,w,t)
    central components = [p,u,v,w]
    upwind blending sensors = [ducrosmod]
    flux function = ausmpw+
}

Note that you have the option to individually centralize primitive variables via central components specification. In the example usage, we have excluded t which stands for temperature. This results in a more stable and smooth scheme compared to full centralization.

CENTRALQ approach results in a low dissipation scheme suitable for LES. However, convective scheme dissipation can be further reduced by using the KEEPBLEND scheme which is the current recommended best practice.

KEEPBLEND Scheme

Instead of blending reconstructed primitive variables at flux faces, KEEPBLEND scheme computes two individual fluxes: i) A non-dissipative kinetic energy and entropy preserving (KEEP) flux (our implementation is approximately entropy preserving), and ii) A dissipative upwind flux. The two fluxes are then blended based on the choice of the upwind sensor. More details about this scheme can be found in Sozer et al.². This scheme will always centralize all flux components, so the central components input is not used. order and limiter only affect the MUSCL reconstruction of the dissipative flux. An example usage is shown below:

convective terms {
    scheme = KEEPBLEND
    # MUSCL reconstruction order [default = 2]
    order   = 2
    limiter = minmod
    upwind blending sensors = [ducrosmod]
    flux function = ausmpw+
}

High Mach Simulations

At high Mach numbers, shocks can develop numerical artifacts, especially if they are not perfectly aligned to the mesh. To remedy this, use the add shock dissipation = true option which employs a shock sensor to detect the shock locations and selectively add extra dissipation.

convective terms {
    # This will employ a shock sensor to add extra dissipation near shocks
    # extra dissipation is computed with Rusanov flux and blended into the
    # flux function selection above
    add shock dissipation = true # [default: false]
}

If performing scale-resolving simulations with ducrosmod or strainmod upwind blending sensors, we recommend setting upwind sensor enstrophy factor = 1 to provide additional robustness at high Mach numbers.

Navier-Stokes Options

This section sets up the option for integrating the Navier-Stokes equations. For explicit time integration, the options are:

navier stokes {
    # Switch between viscous or inviscid
    viscous = true #default
    # Following flow variables are clipped from below and top
    # according to the limits specified. These are not required
    maximum mach        = 2
    minimum density     = 1e-20
    minimum pressure    = 0.1
    maximum pressure    = 1e20
    minimum temperature = 0.1
    maximum temperature = 1e20
}

The minimum and maximum limits of the flow quantities act in various ways. min/max pressure and tempeature are checked during MUSCL reconstruction. If the robust mode is active, these quantities along with maximum Mach number are also used to test for the health of time integration as described in the robust mode section.

For implicit runs, there are a few other inputs, all or most of which can be omitted in favor of using defaults.

navier stokes {
    ...
    # Maximum number of GMRES linear solver iterations
    maximum iterations  = 10 #default
    # Maximum number of Krylov vectors computed with GMRES algorithm.
    subspace size  = 10 # default same as "maximum iterations"
    # Relative (normalized) convergence tolerance for GMRES
    # The linear solver iterations will stop when this level is met
    relative tolerance  = 1.e-3 #default
    # Absolute (non-normalized) converfence tolerance for GMRES
    # The linear solver iterations will stop when this level is met
    absolute tolerance  =  1.e-20 # default
    # If the linear solver residual is below the maximum tolerance
    # adaptive CFL will continue to grow but slowly. If it is exceeded,
    # CFL will be halved
    maximum tolerance = 0.5 # default
    # If the linear solver residual is above the failure tolerance,
    # the implicit step is considered to have failed and the update is rejected
    # The CFL will be set to the minimum specified in the implicit adaptive section
    failure tolerance = 0.8 # default
    # Within each implicit iteration, expected density change in each cell
    # is checked before the linear solver. Cells that will likely exceed the
    # value specified below (20% change allowance) will trigger a reduction
    # in their CFL (hence step size).
    density change threshold = 0.5 # default
    # Similar to above but operates on primary flow variables (p,u,v,w,T) and
    # applied after the linear solver.
    q change threshold       = 0.2 # default
    # The Jacobian calculation algorithm may be chosen separately for Navier-Stokes
    # and turbulence. To override the setting specified in the implicit input section,
    # specify the equation below.
    jacobian method = analytical # analytical (default) or numerical
}

Turbulence Options

LAVA unstructured can currently use Spalart Allmaras (SA) turbulence model and it's various variants for RANS and URANS simulations.

turbulence {
    # Available options: [sa,sa-des,sa-ddes,sa,sa-qcr,sa-qcr-rc,sa-rc]
    model = sa
    # If using qcr variations, pick which version (options: 2000, 2013)
    qcr version = 2000
    # Maximum number of GMRES linear solver iterations
    maximum iterations  = 10 # default
    # Maximum number of Krylov vectors computed with GMRES algorithm.
    subspace size  = 10 # default same as "maximum iterations"
    # Relative (normalized) convergence tolerance for GMRES
    # The linear solver iterations will stop when this level is met
    relative tolerance  = 1.e-3 # default
    # Absolute (non-normalized) convergence tolerance for GMRES
    # The linear solver iterations will stop when this level is met
    absolute tolerance  =  1.e-20 # default
    # If the linear solver residual is below the maximum tolerance
    # adaptive CFL will continue to grow but slowly. If it is exceeded,
    # CFL will be halved
    maximum tolerance = 0.5 # default
    # If the linear solver residual is above the failure tolerance,
    # the implicit step is considered to have failed and the update is rejected
    # The CFL will be set to the minimum specified in the implicit adaptive section
    failure tolerance = 0.8 # default
    # Limit ratio of eddy viscosity to laminar viscosity
    viscosity ratio limit = 1e6 #default
    # Limit change of turbulence variable to 20% in each implicit iteration
    # this guides the adaptive implicit stepping
    change threshold = 0.2 # default
    # Lower limit for SA working variable
    minimum nut = 1e-12 # default
    # Turbulent Prandtl number
    Pr = 0.9 # default
    # Turbulent Schmidt number
    Sc = 0.9 # default
    # The Jacobian calculation algorithm may be chosen separately for Navier-Stokes
    # and turbulence. To override the setting specified in the implicit input section,
    # specify the equation below.
    jacobian method = analytical # analytical (default) or numerical
}

For LES, subgrid scale model is also specified in this section.

turbulence {
    # SGS model (options: none, vreman, smag)
    sgs model = vreman
    # SGS constant for Vreman, default is 0.07
    # but for unstructured solver, best practice is using 0.08
    sgs const = 0.08
}

Wall Function

Wall function is supported for both WMLES and RANS simulations. However, the implicit RANS solver doesn't currently account for wall function Jacobians, leading to decreased stability and lower CFL requirements.

wall function {
    # Pick wall function, options: musker (default), reichardts, blendedlog, rough
    type = blendedlog 
    # If using rough wall model, specify the roughness in meters
    # 0.25e-6m : highly polished
    # 0.50e-6m : polished
    # 3.33e-6m : painted
    roughness = 0.5e-6
    # If activates, this will make the wall function switch between 
    # modeled wall stress and no-slip boundary condition depending on 
    # a wall turbulence sensor. When the sensor indicates laminar, 
    # no-slip will be used
    do transition = false # default
    # Adjusts the threshold where the sensor decides between laminar and turbulent
    # The smaller the number, the earlier the transition
    transition threshold = 0.05 # default
}

Rigid Body Motion

Warning

This feature is a preliminary implementation currently under development. It is only intended for developers and early testers.

This section is used to define the motion for different grid domains, which is essential for simulations involving moving parts, such as those using the sliding mesh method. Each grid file specified in the Grid File and Transformations section is assigned a domain ID based on its order in the list (starting from 1).

You can define multiple motion groups, each identified by an index (e.g., rigidbodymotion_1, rigidbodymotion_2). Each motion group is then assigned to one or more domains.

Two types of rotational motion are supported: constant angular speed and sinusoidal oscillation.

Constant Speed Rotation

Below is an example of setting up a constant speed rotational motion (default):

rigidbodymotion_1 {
    # List of domain IDs this motion applies to
    domains = [1]
    # Type of motion, "rotation" is currently supported
    type = rotation
    # Rotation axis (will be normalized)
    axis = [0,1,0]
    # Center of rotation
    center = [0,0,0]
    # Angular speed in radians per second
    angular speed = 10.0
}

Sinusoidal Rotation

For oscillating motions, you can define a sinusoidal rotation. The angular speed parameter is replaced by amplitude and frequency.

rigidbodymotion_1 {
    # List of domain IDs this motion applies to
    domains = [1]
    # Type of motion
    type = rotation
    # Kind of rotation
    kind = sinusoidal
    # Rotation axis (will be normalized)
    axis = [0,1,0]
    # Center of rotation
    center = [0,0,0]
    # Amplitude of oscillation in radians
    amplitude = 0.1745 # 10 degrees
    # Frequency of oscillation in Hz (revolutions per second)
    frequency = 5.0
}

The parameters are: - domains: A list of domain IDs that will be subjected to this motion. Each domain can only be assigned to one motion group. - type: The type of motion. Currently, only rotation is a valid option. - kind: Kind of motion. For rotation two kinds are supported, constant (default) and sinusoidal - axis: A 3D vector defining the axis of rotation. This vector will be normalized by the solver. - center: A 3D vector specifying the center point of the rotation. - angular speed: The speed of rotation in radians per second (for constant speed rotation). - amplitude: The amplitude of oscillation in radians (for sinusoidal rotation). - frequency: The frequency of oscillation in Hz (for sinusoidal rotation).

These motion definitions are referenced in other sections, like sliding interface, to link a physical motion to a specific part of the simulation domain.

Warning

Currently, only rotational motion (constant and sinusoidal) is supported.

Sliding Mesh

Warning

This feature is a preliminary implementation currently under development. It is only intended for developers and early testers.

Sliding mesh method can be used for efficient simulations of components in rotation to each other, e.g., turbo-machinery applications. Multiple grids corresponding to the domains of different rotation rates (or stationary domain) need to be specified (see grid section).

sliding interface_1 {
    bc pair      = [4,8]
    motion pair  = [1,0]
    run tests    = false # default
    # This will output the ring/panel mesh at the interface in vtk format
    write_interface_mesh = false #default
}

bc pair lists the two boundary condition regions on either side of the interface. motion pair specifies the motion regions corresponding to the bc pair. Here, motion 0 means stationary. See the rigid body motion section for motion specification details.

Note that the boundary conditions for the bc pairs need to be set as sliding (see boundary conditions section).

Initial Conditions

Initial condition specification is optional. If not provided, the code will simply initialize the flow to reference conditions specified in Reference Conditions.

If a different initialization is needed, it can be specified through a numbered list:

# Specify pressure, set the rest using reference conditions
IC_1 {pressure = 5e5}
# Specify temperature, set the rest using reference conditions
IC_2 {temperature = 500}
# If density is specified, another thermodynamic quantity (either pressure or temperature)
# also needs to be specified
IC_3 {density = 2; temperature = 300}
# Velocity vector can be specified
IC_4 {velocity = [100,10,0]}
# Or Mach number along with alpha and beta angles
IC_5 {Mach=0.5; alpha=6; beta=2}
# Specify mass fractions for multi-species
IC_6 {mass fractions = [1,0,0]}

All the above statements will apply the conditions listed in the entire domain. To initialize a sub-volume differently, you can specify a custom box region with position vectors of the two opposing corners:

IC_1 {corner_1=[0,0,0]; corner_2=[1,1,1]; pressure=5e5}

The example above will first apply the reference conditions in the entire domain. Then it will overwrite the pressure in the unit box region specified. The IC's are applied sequentially, i.e. the entries further down in the IC's list will overwrite the previous ones if there is overlap.

A restart file (that is consistent with the current mesh) can also be used as an initial condition via:

IC_1 {
    file = initialCondition.h5
    # Optionally, you can apply smoothing to the file being read as initial condition via:
    apply smoothing passes = 2 # default is 0
}

This feature is useful when starting a simulation from another solution interpolated onto the current grid (which can be done using unst_interpolator tool)

Warning

Note that if the file method is used, no other initial conditions can be specified.

Artificial Boundary Layer

Artificial boundary layer options will linearly scale the velocity down to zero at walls, creating a less harsh start-up. For isothermal walls, same linear scaling can be applied to the temperature field.

Warning

This option, while available, should be avoided for large scale LES cases, as it requires distance function to be computed for the entire volume which is otherwise not needed.

artificial boundary layer {
    apply velocity = true
    apply thermal  = true
    # Distance from the wall (in meters) that the linear scaling is applied to
    height = 1
    wall temperature = 500 # K
}

Sponge

Sponge zones can be used to dampen the flow in specified rectangular regions. In these zones, the order of accuracy will be reduced.

sponge_1 {
    corner1 = [-10.0, -10.0, -10.0]
    corner2 = [-5.0, 10.0, 10.0]
}

Multiple sponge regions can be defined by incrementing the number after sponge_

Time Averaging

This section controls whether time averaging of the flow field is requested (if the section is present) and sets the starting point.

time averaging {
    # Start time averaging at a certain step, physical time, or CTU
    # (Convective Time Unit, i.e. non-dimensional time = physical time / reference time).
    # Only one of the following three options can be specified.
    # If none are specified, averaging begins from the start of the simulation.

    #start step = 1000
    #start time = 1.0
    start ctu  = 50

    # Activate time averaging of surface quantities (default: true)
    do surface = true
    # Activate time averaging of volume quantities (default: true)
    do volume  = true

    # Sampling interval: take an averaging sample only every N time steps.
    # Default is 1 (sample every step).
    sampling interval = 1
}

Exponential Moving Average

By default, LAVA uses a cumulative average (equal weight to all past samples). For long unsteady runs it can be useful to instead use an exponential moving average (EMA), which gives a decaying weight to older data so that the average tracks recent flow statistics.

time averaging {
    start ctu = 50

    # Switch from cumulative average to exponential moving average
    do exponential moving average = true

    # Window size W (in CTU). Data older than W has approximately 3% residual weight.
    # Specify in CTU (non-dimensional, CTU = physical time / reference time):
    exponential moving average window ctu = 10

    # Alternatively, specify the window directly in physical time units (seconds):
    # exponential moving average window = 0.05
}

The EMA blending weight at each sample step is:

\[ \alpha = \min\!\left(1,\; \frac{c}{W} \cdot \Delta t_\text{elapsed}\right), \qquad c = 3.5 \]

where \(W\) is the window size in physical time and \(\Delta t_\text{elapsed}\) is the elapsed physical time since the last averaging sample. At each step the average is updated as:

\[ \overline{q} \leftarrow \alpha\, q + (1 - \alpha)\, \overline{q} \]

In the continuous limit, data from time \(\tau\) ago carries weight \(\propto e^{-c\,\tau/W}\). At exactly one window back (\(\tau = W\)) the residual weight is \(e^{-3.5} \approx 3\%\).

Trip Lines

Trip lines exert a penalty method source term that attempts to trigger boundary layer transition to turbulence. The source terms are applied only at select cells within certain distance of the lines. Lines are also broken into dash/dot or on/off pattern that apply the trip terms in an intermittent fashion along the line, similar to how trip dots are commonly implemented in reality.

# This determines the strength of the trip lines
# The lower the value, the stronger the trip source penalties are
# In the lower limit, a value of 1 would bring the flow to full stop
# where trip is active
trip time relaxation = 40. # a factor of dt
trip_line_1 {
    # dxscale should indicate the local grid resolution (where this trip line is)
    # with respect to the dx value set at the reference section
    dxscale = 1
    # Begin/end coordinates of the line segment
    begin = [0,1,1]
    end   = [1,2,3]
    # The segment variables below are relative to the local grid spacing
    # In this example, the trips will be on for approximately 2 cells
    # and off for approximately 3 cells along the line
    # The on segments will mark cells whose centers are within 0.9 * dx
    # from the line
    # These variables can be skipped in favor of defaults
    segment length = 2 # default
    segment skip   = 3 # default
    segment radius = 0.9 # default
    # The straight trip line specified here will be projected to the
    # surface using projection method. By default, the line will be projected
    # to the closets surface
    # Options are : none, closest, vector
    projection method = closest # default
    # max distance is an optimization parameters, entering the something 
    # slightly larger than the expected farthest distance (from trip line to the geometry) 
    # will improve start-up performance
    projection max distance = 3
}
trip_line_2 {
    dxscale = 4
    begin   = [1,0,0]
    end     = [2,0,0]
    # Another projection method is vector, for which the line 
    # is projected along a specified vector direction until a surface is
    # encountered
    projection method = vector
    projection vector = [0,1,0] # y-axis
}
# Circular trip curve
trip_line_3 {
    dxscale = 2
    # Instead of a straight line, the trip location can also be drawn as a circle
    # with the specified center coordinates and radius
    center  = [1,0,0]
    radius  = 2
    # The circle can also be projected to the surface
    projection method = closest 
    projection max distance = 1
}

Tip

In typical practice, it is recommended to use defaults for segmenting parameters. In addition, the inputs may be made more compact by inlining each trip, for example:

trip time relaxation = 40. 
trip_line_1 {dxscale=1;begin=[0,1,1];end=[1,2,3];projectionmaxdistance=3}
trip_line_2 {dxscale=4;begin=[1,0,0];end=[2,0,0];projectionmethod=vector;projectionvector=[0,1,0]}
trip_line_3 {dxscale=2;center=[1,0,0];radius=2;projectionmethod=closest;projectionmaxdistance=1}

Tip

The solver will output the final projected trip lines (in tecplot and vtk line formats) into the trips directory under the run directory. It is highly recommended to verify the trip setup by plotting the lines overlaid on the surface geometry.

Output

All the output settings are organized under output {} section.

Flow Stats

Controls interval of sampling flow statistics such as integrated kinetic energy, maximum Mach number and integrated domain volume (see Monitor section for descriptions).

output {
    flow stats {
        # Defaults to monitor interval by default
        interval = 10
    }
}

Monitor

Monitor output tracks typically high frequency information such as time step, CFL, residuals, integrated aerodynamics loads and performance metrics among others.

If monitor section is not included in the input file, you will get a default set of variables output to screen every time steps (or every pseudo step for implicit runs). The variable list given here only affects the screen output while every other available variable is always output to the .iter file.

The .iter file can be plotted using the tool plotHist.py

output {
    monitor {
        # Output interval applies both to screen and file (.iter) file outputs.
        interval = 10
        # Variable list only affects the screen output
        variables = [nt,time,ctu,dt,maxcfl,maxmach,cl_airplane,wallt,mups,elapsed]
        # When restarting, the code looks for a copy of the .iter file in the restart directory
        # corresponding to the time step you are restarting from. If the .iter file is missing,
        # by default the code will quit with an error message. The flag below, if activated, will
        # instead start a new .iter file from scratch
        reset if missing = false # default
        # The monitor keeps track of some averaged performance metrics. This entry will set the interval
        # for averaging. By default, it is set to 10xinterval (specified above)
        avg interval = 100
    }
}

List of available monitor variables:

Variable	Description	Units
nt	Time step
ns	Pseudo time step
time	Physical simulation time	[s]
ctu	Convective time unit	[non-dimensional time]
dt	Time step size	[s]
maxcfl	Maximum CFL number in the domain
maxmach	Maximum Mach number in the domain
integratedmass	Integrated mass of fluid in the domain	[kg]
integratedvolume	Integrated volume of the domain	[m^3]
integratedke	Integrated kinetic energy in the domain	[J]
nondim_integratedke	Non-dimensional integrated kinetic energy
integratedxmomentum	Integrated x-momentum	[kg*m/s]
integratedymomentum	Integrated y-momentum	[kg*m/s]
integratedzmomentum	Integrated z-momentum	[kg*m/s]
integratedetotal	Integrated total energy in the domain	[ (kg/m^3)(m^2/s^2) ]
nondim_integratedetotal	Non-dimensional integrated total energy
minpressure	min pressure in the domain	[Pa]
maxpressure	max pressure in the domain	[Pa]
minuvelo	min u-velocity in the domain	[m/s]
maxuvelo	max u-velocity in the domain	[m/s]
minvvelo	min v-velocity in the domain	[m/s]
maxvvelo	max v-velocity in the domain	[m/s]
minwvelo	min w-velocity in the domain	[m/s]
maxwvelo	max w-velocity in the domain	[m/s]
mintemperature	min temperature in the domain	[K]
maxtemperature	max temperature in the domain	[K]
mups	Average MUPS, averaged through the log interval
mups_min	Minimum MUPS through the monitor log interval
mups_max	Maximum MUPS through the monitor log interval
mups_long	Average MUPS, averaged through the avg interval
mups_long_min	Minimum MUPS through the monitor avg interval
mups_long_max	Maximum MUPS through the monitor avg interval
wallt	Wall time per time step, averaged through the log interval	[s]
wallt_min	Minimum wall time through the monitor log interval	[s]
wallt_max	Maximum wall time through the monitor log interval	[s]
wallt_long	Wall time (s), averaged through the avg interval	[s]
wallt_long_min	Minimum wall time through the monitor avg interval	[s]
wallt_long_max	Maximum wall time through the monitor avg interval	[s]
elapsed	Total elapsed wall time from beginning, including the restart history	[hr]
memory	CPU host memory usage on root MPI rank	[GB]
memory_gpu	GPU device memory usage on root MPI rank	[GB]

Note

MUPS: Millions of Updates per Second, i.e. millions of cells advanced through a full time step per second

Integratrated loads, as described in the Loads section, are also output via the monitor. Each loads integration group, corresponding to the its specied name will produce the following monitor output variables:

Variable	Description
cx_name	Total x-force coefficient
cy_name	Total y-force coefficient
cz_name	Total z-force coefficient
cx_inv_name	Inviscid x-force coefficient
cy_inv_name	Inviscid y-force coefficient
cz_inv_name	Inviscid z-force coefficient
cx_visc_name	Viscous x-force coefficient
cy_visc_name	Viscous y-force coefficient
cz_visc_name	Viscous z-force coefficient
cmx_name	Total x-moment coefficient
cmy_name	Total y-moment coefficient
cmz_name	Total z-moment coefficient
cmx_inv_name	Inviscid x-moment coefficient
cmy_inv_name	Inviscid y-moment coefficient
cmz_inv_name	Inviscid z-moment coefficient
cmx_visc_name	Viscous x-moment coefficient
cmy_visc_name	Viscous y-moment coefficient
cmz_visc_name	Viscous z-moment coefficient
cl_name	Total lift coefficient
cd_name	Total drag coefficient
cs_name	Total side force coefficient
cl_inv_name	Inviscid lift coefficient
cd_inv_name	Inviscid drag coefficient
cs_inv_name	Inviscid side force coefficient
cl_visc_name	Viscous lift coefficient
cd_visc_name	Viscous drag coefficient
cs_visc_name	Viscous side force coefficient

If running the implicit solver, there are a few additional monitor variables available:

Note

For residual norms ICOMP index below refers to the equation number: 0: Continuity, 1-2-3: Momentum, 4: Energy, >=5: Species mass fractions For DQ norms, ICOMP index below refers to the primitive variables: 0: Pressure, 1-2-3: U-V-W velocity components, 4: Temperature, >=5: Species mass fractions

Variable	Description
l0resICOMP	L0 Residual norm (absolute)
l2resICOMP	L2 Residual norm (absolute)
l0resrICOMP	L0 Residual norm (normalized by first step value)
l2resrICOMP	L2 Residual norm (normalized by first step value)
l0dqICOMP	L0 Primitive variable change norm (absolute)
l2dqICOMP	L2 Primitive variable change norm (absolute)
l0dqrICOMP	L0 Primitive variable change norm (normalized by first step value)
l2dqrICOMP	L2 Primitive variable change norm (normalized by first step value)
l0turbres	L0 Residual norm for turbulence equation (absolute)
l2turbres	L2 Residual norm for turbulence equation (absolute)
l0turbresr	L0 Residual norm for turbulence equation (normalized by first step value)
l2turbresr	L2 Residual norm for turbulence equation (normalized by first step value)
lambdap	Under-relaxation factor for Navier-Stokes equations
lambdatur	Under-relaxation factor for turbulence equation
lambdals	Under-relaxation factor for line search
linits	Number of linear solver (GMRES) iterations for Navier-Stokes equations
linres	Linear solver (GMRES) residual for Navier-Stokes equations
linitsturb	Number of linear solver (GMRES) iterations for turbulence equation
linresturb	Linear solver (GMRES) residual for turbulence equation

Loads

In this section, you can define several loads integration groups which can integrate forces and moments over specified components of the geometry.

output {
    loads {
        # Interval at which the loads are computed, this is common for all integration groups
        interval      = 10
        # Following 3 inputs apply to all the groups by default
        # But you can optionally specify them different for any individual group
        # Center around which moments are computed
        moment center = [33.,0.0,4.5]
        # Provide non-dimensional force and moments
        normalize     = true
        # Use gauge pressure (p-p_ref) in pressure forces
        gauge         = true
        # If true, this flag will output each surface component integrations separately
        # (in addition to the groups listed below)
        do individual comps = true
        group_1 {
            # Group name
            name          = airplane
            # Which BC zones to include
            include       = [1-32]
        }
        group_2 {
            name = wing
            include = [7,8,10-13]
            normalize     = false
            # Use gauge pressure (p-p_ref) in pressure forces
            gauge         = false
        }
        group_3 {
            name = tail
            include = [14,16]
            moment center = [48.,0.0,4.5]
        }
        ...
    }
}

The code will always output non-normalized force/moment values. If normalize = true is specified, the code will additionally output normalized force/moment coefficients. Reference section values are used for normalization.

The result of the load integrations are output to screen (if requested) and to the monitor file (.iter). Detailed description of the monitor and the loads output variables are provided in the monitor section.

If you only have single loads integration group, you can skip the group_1 subsections:

output {
    loads {
        interval      = 10
        moment center = [33.,0.0,4.5]
        normalize     = true
        gauge         = true
        do individual comps = true
        name          = airplane
        include       = [1-32]
    }       
}

Restart

Restart files, also referred to as checkpoint files, allow you to restart the simulation seamlessly from a given point. You can set the interval of restart file writes as:

output {
    restart {
        # Specify restart file output interval (time step for unsteady, pseudo-step for steady)
        interval = 1000
    }
}

Volume

Volume output produces a CGNS file with the selected flow variables. This will create one mesh link file which contains the grid connectivity information and a solution file for each output timestep.

output {
    volume {
        # Whether to write out an volume file at the start of simulations (initial conditions)
        write initial = false
        # Currently only cgns format is supported
        format        = cgns # default
        # Interval of output (unsteady: time step, steady: pseudo-step)
        interval      = 10000
        # List of variables to output
        variables     = [q,qstats,density,mach,problem,sponge,upwind,cellmask,trip]
    }
}

In addition to the full volume output, you can also output a subset of the domain. Currently, two different types of subset are supported. They are described below along with the supporting options required:

• box: Output all the cells inside a box defined by box vertices 'lo' and 'hi'. The 'lo' parameter defines the minimum x,y,z coordinates of the box, while 'hi' defines the maximum x,y,z coordinates.
• plane: Output all the cells that are intersected by a plane defined by 'normal' and 'origin'. The 'normal' parameter defines the direction perpendicular to the plane (does not need to be normalized), and 'origin' specifies any point that lies on the plane. The 'clip lo' and 'clip hi' parameters can optionally be used to clip the plane output to a specific bounding box region.

  volume {
        # Whether to write out an volume file at the start of simulations (initial conditions)
        write initial = false
        # Currently only cgns format is supported
        format        = cgns # default
        # Interval of output (unsteady: time step, steady: pseudo-step)
        interval      = 10000
        # List of variables to output
        variables     = [q,qstats,density,mach,problem,sponge,upwind,cellmask,trip]

        subset_1 {
            type = box
            lo = [0.1, 0, 1.0]    # Minimum coordinates [x_min, y_min, z_min]
            hi = [1.0, 0.5, 0.1]  # Maximum coordinates [x_max, y_max, z_max]
            interval = 500        # Output interval specific to this subset
        }
        subset_2 {
            type = plane
            normal = [0, 0.5, 0.5]      # Direction vector perpendicular to the plane
            origin = [1.0, 0., 0.1]     # A point located on the plane
            cliplo = [-2.0, -0.5, -1.0] # (Optional) For clipping plane: Minimum coordinates [x_min, y_min, z_min]
            cliphi = [2.0, 0.5, 1.0]    # (Optional) For clipping plane: Maximum coordinates [x_max, y_max, z_max]
            interval = 100              # Output interval specific to this subset
        }
    }

The variables, write initial and format inputs are shared across all subsets. Each subset can have its own output interval that differs from the main volume inteval. Output files will follow the naming convention of sub[subset_id]_[solver].vol.[time_step].cgns.

Tip

For large simulations, using subsets can significantly reduce storage requirements and post-processing time.

Full list of volume output variable options:

Variable	Description
pressure
velocity	All velocity components
temperature
massfractions	Species mass fractions, only supported for multi-species
q	All of the above primitive flow variables
qstats	Time average and statistics of the flow field
density
viscosity	laminar dynamic viscosity
nut	SA turbulence model working variable
eddyviscosity (or mut)
mach	Mach number
distancefunction	Distance from walls
limiter	MUSCL slope limiter for each flow variable
upwind	Upwind/central blending sensor field, 0: central, 1: upwind
volume	Volume of each cell
centroid	All 3 centroid coordinates for each cell
lsmooth	Smoothed length scale for each cell
residuals	All residuals (implicit only)
massresidual	Continuity equation residual implicit only
xmomresidual	X-momentum residual (implicit only)
ymomresidual	Y-momentum residual (implicit only)
zmomresidual	Z-momentum residual (implicit only)
energyresidual	Energy equation residual (implicit only)
nutresidual	SA turbulence model residual
dtau	Pseudo-time step size
sponge	Binary field indicating where sponging is active
trip	Binary field indicating where trips are active
thin	Binary field indicating whether a cell is marked as thin
problem	Binary field indicating whether a cell is marked as problem
problem_gg	Binary field indicating whether a cell failed Green-Gauss check
id	partition id for each cell
tiles	tile id for each cell
domain	domain id for each cell
motion	motion id for each cell (-1: stationary)
cellmask	Cell bit mask

Surface

Surface output produces a CGNS file with the selected flow variables for all boundary condition zones. By default, output uses solution linking: a single mesh file (Unst.surf.meshLink.cgns) is written once at startup and each solution file links to it rather than duplicating the geometry data. This significantly reduces file sizes for long runs.

output {
    surface {
        # Whether to write out an initial file at the start of simulations (initial conditions)
        write initial = false
        # Interval of output (unsteady: time step, steady: pseudo-step)
        interval      = 10000
        # List of variables to output
        variables     = [q,qstats,density,mach,cf,yplus]
        # Write geometry and solution in a single self-contained file instead of linking.
        # Default: false. Set to true to output the solution and a linked mesh file separately.
        split grid and solution = false
    }
}

BC Statistics Output

At every surface output step, area-averaged boundary condition statistics are appended to surface/bc_stats.iter, a continuous binary file in the same format as monitor/unst.data.iter. A human-readable mapping of BC indices to names is written once to surface/bc_name_map.txt.

The following quantities are recorded per boundary condition (prefixed BC1_, BC2_, etc.):

Column	Description
step, time	Simulation step and time
BCN_area	Total boundary face area (m²)
BCN_pressure	Area-averaged pressure (Pa)
BCN_u, BCN_v, BCN_w	Area-averaged velocity components (m/s)
BCN_temp	Area-averaged temperature (K)
BCN_y0, BCN_y1, ...	Area-averaged species mass fractions (multi-species only)
BCN_mach	Area-averaged Mach number
BCN_mdot	Area-averaged mass flux (kg/(m²·s))
BCN_tauw_x/y/z	Area-averaged wall shear stress components (Pa)
BCN_heatflux	Area-averaged heat flux (W/m²)
BCN_hflux_min, BCN_hflux_max	Min/max heat flux over the boundary (W/m²)

When mach-driven pressure BCs are active, per-BC pressure history is written to surface/mach_driven_pressure_BCN.iter with columns step, mach_ave, pressure_ave.

Both bc_stats.iter and mach_driven_pressure_BCN.iter support restart continuation: history is automatically preserved across restart segments.

Full list of surface output variable options:

Variable	Description
pressure
velocity	All velocity components
temperature
massfractions	Species mass fractions, only supported for multi-species
q	All of the above primitive flow variables
density
mach	Mach number
viscosity	laminar dynamic viscosity
conductivity	thermal conductivity
specificheat	Cp
dtau	Pseudo-time step size
mdot	Mass flux (kg/(m^2.s), positive: inflow)
qdot	Heat flux (W/m^2, positive: inflow)
tauw	Wall shear stress (vector components)
cf	Surface skin friction coefficient (vector components)
pcoef	Surface pressure coefficient
yplus	Wall yplus value
nut	SA turbulence model working variable
facemask	Face bit mask
cellmask	Cell bit mask
qavg	Time-averaged primitive flow variables
qstats	Time average and statistics of the flow field
tauwavg	Time average of wall shear stress (vector components)
yplusavg	Time average of wall yplus value
turbulence	Turbulence sensor

Iso-Surfaces

In-situ extraction of iso-surfaces for post-processing and visualization is supported. Iso-surfaces are defined by specifying a contour value for a chosen field, and the solver will output the corresponding isosurface as a triangulation at a user-defined interval.

isosurfaces_1 {
    isovalue     = 0.5
    contour      = vortmag
    colorfield   = umag
    interval     = 100
    format       = vtk
    clipbox_1 {
        lo     = [0.0, 0.0, 0.0]
        hi     = [1.0, 1.0, 1.0]
        invert = false
    }
}
isosurfaces_2 {
    contour      = cutplane
    normal       = [0, 1, 0]         # Normal vector of the plane (e.g., y-normal)
    point        = [0.0, 0.5, 0.0]   # A point on the plane
    colorfield   = densitygradmag
    interval     = 100
}
isosurfaces_3 {
    contour     = cutplane
    point       = [9.75, 0, 0]
    normal      = [1, 0, 0]
    format      = vtk
    interval    = 100
    colorfield  = fullq
    derivedvars = [mach, density]
}
isosurfaces_4 {
    isovalue        = 300
    contour         = q
    qcompcontouring = 4 # temperature component of the q vector for 3D simulations
    colorfield      = q
    qcompcoloring   = 0 # pressure component of the q vector
    interval        = 100
    format          = vtk
}

Parameters

• isovalue: The value of the chosen field at which to extract the iso-surface (e.g., isovalue = 0.5).
• contour: The field to contour. Options include:
  • vortmag (vorticity magnitude, default)
  • qcriterion
  • mach
  • densitygradmag
  • q (a specific component of the q vector. The option qcompcontouring is then required to specify the component index)
  • cutplane (for planar cuts; requires normal and point)
• colorfield: Field used to color the iso-surface. Options include:
  • umag (color with velocity magnitude, default)
  • vortmag
  • qcriterion
  • errormetric
  • densitygradmag
  • q (color with a specific component of the q vector. The option qcompcoloring is then required to specify the component index. If the option qcompcoloring is missing, then fullq is used instead)
  • fullq (color with all components of the q vector)
• interval: Output interval in time steps (e.g., interval = 100).
• format: Output file format. Buffered Options: vtk (default), vtp, cgns, fwh. Non-buffered options: xdmf, vtk, vtp.
  • vtp writes VTK XML PolyData (.vtp) files with raw binary appended data.
• cliplo, cliphi: (Optional, deprecated) Lower and upper bounds for spatial clipping of the iso-surface, specified as coordinate arrays. Defaults to the full domain. Equivalent to a single clipbox_1 { lo=...; hi=...; invert=false }; emits a deprecation warning at run time. Cannot be combined with indexed clipbox_N sub-sections. Prefer clipbox_1 in new inputs.
• clipbox_N: (Optional) Indexed axis-aligned box clip(s). Up to 4 boxes are supported per iso-surface (clipbox_1, clipbox_2, ...). Each box has lo, hi, and an optional invert flag (default false). When invert=false the box keeps points inside it; when invert=true it keeps points outside it. All configured boxes are AND'd together, and AND'd with all clipfrustum_N regions. Example — keep points inside a large box but excluded from a small interior box:

  clipbox_1 { lo = [0.5, 0.0, 0.3];  hi = [1.0, 0.5, 0.7];  invert = false }
  clipbox_2 { lo = [0.7, 0.2, 0.45]; hi = [0.8, 0.3, 0.55]; invert = true  }
  

• clipfrustum_N: (Optional) Indexed frustum (or hollow frustum) clip(s). Up to 4 frustums are supported per iso-surface (clipfrustum_1, clipfrustum_2, ...). The frustum axis runs from lo to hi; radius is the cross-section radius at lo, and the optional radiushi is the radius at hi (defaults to radius, i.e., a cylinder). The optional radiusinner and radiushiinner carve out a coaxial hollow region (defaults: 0). An optional invert flag (default false) flips the keep/discard sense for that frustum: with invert=true points inside the frustum are discarded. All frustums are AND'd together and AND'd with any box clips. Example — a hollow cylinder along the y-axis between y=0.1 and y=0.4 with outer radius 0.18 and inner radius 0.08:

  isosurfaces_4 {
      isovalue   = 10
      contour    = vortmag
      colorfield = umag
      interval   = 100
      clipfrustum_1 {
          lo            = [0.75, 0.10, 0.50]
          hi            = [0.75, 0.40, 0.50]
          radius        = 0.18
          radiusinner   = 0.08
          # radiushi      = 0.18     # optional, defaults to radius (cylinder)
          # radiushiinner = 0.08     # optional, defaults to radiusinner
          # invert        = false    # optional, default false
      }
  }
  

• normal, point: (Optional) For planar cuts, specify the normal vector and a point on the plane.
• qcompcontouring: Required for contour = q, selects the q component to contour.
• qcompcoloring: Required for colorfield = q, selects the q component for coloring.
• derivedvars : (Optional) For colorfield = fullq or colorfield=q without any qcompcoloring specification, the output state vector can be appended with derived variables. Available options are any combinations of:
  • mach
  • density

Multiple Iso-Surfaces

You can define multiple iso-surfaces by incrementing the section name, e.g., isosurfaces_1, isosurfaces_2, etc.

Notes

• If contour = cutplane, you must provide both normal and point.
• If format = xdmf, parallel output is enabled by default (only recommended for very large iso files).
• If colorfield = q and qcompcoloring is not specified, the full Q vector is used for coloring.
• There is a parallel/serial option exposed to the user, but we recommend to not set that and allow the code to determine the best setting.
• Buffered output for isosurfaces is controlled by the outputbufferintervalisos = 5000 variable in the output section. No need to set that variable as the code will automatically determine a best buffering interval.
• For qcompcontouring and qcompcoloring, the component of the state vector follows this 0-index order: [p u v (w) T (Y1 Y2 Y3 ...)]. Here, note that the field w only exists when running the 3D executable.

Note

LAVA unstructured currently only supports cgns output format for isosurfaces. In addition, the errormetric contour or color fields are not supported.

Warning

Isosurface extraction is currently only supported on GPU runs.

Surface Probe

Surface probes sample the flow field on to user supplied triangulated surfaces:

Note

LAVA unstructured currently only supports cgns and fwh output format for surface probes.

output {
    surface probe_2 {
        # Specify file, the format is determined by extension
        # Supported formats: tri, fro, cgns
        file = probe_triangulation.cgns
        interval = 100 # Interval of output
        # Where to sample data, either face center or nodes (vertices)
        data location = face # default, options: face, node
        # Optionally transform the probe triangulation. The operations are applied in the order: 1. rotation, 2. scale, 3. translate
        rotate = [90.,0.,0.]
        scale = [2.,1.,1.]
        translate = [1., 1. ,1.]
        # Select output format, available options: vtk, cgns
        format = cgns
        output cell status = true # Output cell status information like limiters, sponge, upwind
        output integrated  = true # Output integrated quantities like fluxes through the surface probe
    }
}

By default, surface probes will sample primitive flow variables (p,u,v,w,T,Y). In addition, an inout variable will be written (whether the sampling point was able to find a donor cell or not). Cell status information like problem cells, limiter, upwind sensor and sponge can be requested by setting the option output cell status to true.

Integrated quantities like mass flux can be requested by setting the option output integrated to true.

Averaged quantities can be requested by setting the option output qstats to true.

output {
    surface probe_2 {
        ...
        output qstats = true # Output averaged quantities sampled on to the probe
    }
}

The results are written out in two parts, a grid link file that is created during the job startup and solution files that are written out at the requested sampling interval.

For high frequency sampling of acoustic surfaces, a specialized .fwh format is implemented that only writes the sampled data without mesh information.

output {
    surface_probe_1 {
        file = FWH_D0.67_35.i.tri
        format = fwh
        interval = 10
    }
}

Line Probe

lineprobe_1 {
    start = [-1.000, 0.0,0];
    end = [1.000,0.0,0];
    number of points = 512;
    interval = 20;
    }

Warning

By default, each line probe output is written to an ASCII formatted file with one file per time step. This can quickly generate numerous files, potentially causing file system issues.

To mitigate this, a binary file format can be used by setting write binary to true. This creates a single binary file per line probe that contains data for all timesteps, significantly reducing the number of files created during your simulation.

output {
    line probe_1 {
        ...
        write binary = true # Output line probes to a binary formatted file
    }
}

For visualization and analysis, you can use the plotLineProbe.py tool, which is specifically designed to read and process these binary files.

When restarting simulations, the code handles existing binary line probe files:

1. It checks if the current line probe definition matches the existing file
2. If definitions match and restart time step is consistent, it simply appends new data
3. If definitions match but restart time step differs, it creates a backup of the original file and prepends existing data up to the restart point in a new file
4. If definitions don't match, it creates a copy of the existing file and starts a new file with the current definition

This robust handling ensures data integrity across simulation restarts while preserving your previous results.

Averaged quantities can be requested by setting the option output qstats to true.

output {
    line probe_1 {
        ...
        output qstats = true # Output averaged quantities sampled on to the probe
    }
}

Wall-Normal BL Probe

Line probes support a wall normal mode for boundary layer profiling. Instead of specifying a start/end pair, you supply one or more seed coordinates; the solver projects each seed onto the nearest no-slip or wall-function boundary face and grows a wall-normal ray of sample points outward to a specified height.

Seed specification — a single seed can be given inline with start; multiple seeds are read from a plain ASCII file (one x y z per line):

output {
    line probe_2 {
        wall normal = true
        file        = seeds.dat   # one "x y z" per line; omit to use start = [x,y,z]
        height      = 0.1         # outer extent of each ray in solver length units
        interval    = 500
        output qstats = true
        name          = bl_probe
    }
}

Ray discretization — choose one of two modes:

Uniform — evenly spaced points including the foot at h = 0:

        number of points = 50

Geometric — first cell height ds0, each subsequent step multiplied by growth rate (default 1.2) until height is reached; the final point is snapped to height. The actual number of points is determined automatically.

        initial spacing = 5e-6
        growth rate     = 1.2     # optional, default 1.2

Optional keys:

Key	Type	Default	Description
wall bc names	string list	all wall BCs	Restrict BVH projection to BCs matching these names
flip normal	bool	false	Flip the wall-outward normal direction
output qstats	bool	false	Include time-averaged quantities (requires volume time averaging)
name	string	line_probeN	Base name for the output file

Output file format:

BL probe files are written in Tecplot ASCII multi-zone format. One zone is created per seed; the zone title carries the 0-based seed index (T="seed_0", T="seed_1", …). Each zone header stores the projected foot coordinates as AUXDATA:

AUXDATA key	Description
x_foot, y_foot, z_foot	Projected wall face centroid coordinates

Data columns (one row per sample point):

nt  time  x  y  z  [flow vars]  [qstats vars]  inout  height  tau_w_inst  rho_w_inst  mu_w_inst  [tau_w_avg  rho_w_avg  mu_w_avg]

• height — wall-normal distance from the foot (0 at the projected wall face centroid)
• tau_w_inst / tau_w_avg — instantaneous and time-averaged wall shear stress magnitude at the foot (constant within a zone)
• rho_w_inst / rho_w_avg — density at the foot, instantaneous and time-averaged
• mu_w_inst / mu_w_avg — dynamic viscosity at the foot, instantaneous and time-averaged
• Avg columns (tau_w_avg, rho_w_avg, mu_w_avg) are present only when output qstats = true

On the first output step, two Tecplot equation files are written to the line_probe/ directory:

• bl_inner_scale_inst.eqn — computes u_tau, nu_w, y+, u+, v+, w+ from the instantaneous foot columns
• bl_inner_scale_avg.eqn — same using _avg columns; also adds <uu>+ through <vw>+ from qstats

Load via Data → Alter → Specify Equations → Load in Tecplot.

Point Probe

Point probes will sample the flow state at a number of given points. A group of points can be specified at once using an ascii file that lists the x y z coordinates of each point in order (space delimited):

output {
    point probe_1 {
        file = ascii_file_with_xyz_coords.dat
        interval = 1 # Interval of sampling
        scale = 1 # Option to scale the coordinates
    }
}

Results will be written to an ascii file in the point_output directory, where each probe section and each variable will get its own file.

Alternatively, points can be output to the iter file and be plotted using plotHist.py tool by:

output {
    point probe_1 {
        ...
        write to monitor = true # write the probes to the .iter file
    }
}

Point probes will sample primitive flow variables (p,u,v,w,T,Y_i) along with an inout variable indicating whether the query point was detected inside the fluid domain or out. Optionally, you can output the Mach number as:

output {
    point probe_1 {
        ...
        do mach = 1 # Sample Mach number as well
    }
}

Averaged quantities can be requested by setting the option output qstats to true.

output {
    point probe_1 {
        ...
        output qstats = true # Output averaged quantities sampled on to the probe
    }
}

Boundary Conditions

Boundary conditions are specified in the form of a numbered list, where the order must match the boundary zones in the CGNS grid file:

BC_1 {type=farfield}
BC_2 {type=wall}
BC_3 {type=symmetry}

The order of the boundary conditions (BC's) needs to be consistent with the order of BC's zones in the grid file. A list of BC zones in the grid file can be obtained using the tool unst_gridInfo.

Available BC types are wall, wallfunction, symmetry, centerline, inlet, outlet, farfield, exact, nozzle, distributed and sliding. Each is explained below with their options.

Wall

Wall type boundaries can either impose no-slip or slip conditions, specified with the kind statement:

The wall BC is treated as adiabatic by default. You can change the thermal behavior by either specifying a temperature (for isothermal wall) or a heat transfer rate (for fixed heat transfer wall). The heat transfer rate qdot should be given in [W/m^2] units. A positive qdot indicates energy being added to the domain (heating).

# Isothermal wall fixed at 300K
BC_1 {type=wall; temperature=300}

# Constant heat transfer wall with 1000W/m^2
BC_2 {type=wall; qdot=1000}

Conjugate fluid/solid heat conduction (solved in 1D rays in the wall normal direction, dicretized through the solid thickness) modeling can be activated through wall boundary conditions. This model needs wall material properties and thickness as inputs.

BC_1 {
    type = wall
    temperature = 300 # Used as initial condition
    conjugate solid = Steel_ASTM_A572_Grade50.sol
    # Number of points (uniformly distributed) in the solid ray
    solid number of points = 300
    solid thickness = 0.02
    solid output stride = 10 # Output every 10th conjugate ray
}

The format of the conjugate solid material file is explained here.

Wallfunction

Wall function is an alternative wall boundary condition that applies a wall model:

BC_1 {type=wallfunction}

The default behavior of the wall function is controlled through the wallfunction section in the input file. However you can customize the wall function type and roughness settings for each BC as:

BC_1 {type=wallfunction; wf_type=rough; roughness=1.e5}

Symmetry

Symmetry BC has no additional inputs.

BC_1 {type=symmetry}

Centerline

For axisymmetric simulations:

BC_1 {type=centerline}

Please see grid transformation section for more details on how to set up an axisymmetric simulation.

Farfield

This is characteristics based inlet/outlet farfield condition for all speeds. If no options are given, it will use the refence conditions as farfield values.

BC_1 {type = farfield}

You can also set the farfield conditions to be different than reference conditions:

BC_1 {
    type = farfield
    pressure = 1e5
    temperature = 300
    # density = 1.2
    velocity = [100,0,0]
}

You can specify either two of the thermodynamic quantities: pressure, temperature and density.

Inlet

Inlet is a flexible type that will activate the proper boundary condition based on what is specified. Several different inlet implementations are available as detailed below. Before specific details, here are some common inlet options that are shared:

BC_1 {
    type=inlet;

    # Mass fractions are only needed for multi-species simulations.
    massfractions=[0.233, 0.767]  # For 2 species

    # For turbulence models
    eddyviscosityratio=3.0;

    # Optional sponge layers to dampen reflections
    spongelayers=2;


    }

Velocity inlet

If velocity is specified, velocity inlet BC is activated:

BC_1 {
    type=inlet;
    # Required
    velocity=[100,0,0];   # m/s

    # One of these options is required for subsonic inlet
    density=1.225;        # kg/m^3
    #pressure=101325;      # Pascal
    #temperature=300;      # Kelvin
}

Stagnation inlet

Stagnation conditions (total pressure/temperature) can be used as:

BC_1 {
    type=inlet;
    totalpressure=101325;     # Pascal
    totaltemperature=300;     # Kelvin
}

Mass flow inlet

BC_1 {
    type=inlet;
    mdot=0.1;             # kg/(m^2.s) (positive for inflow)

    # One of these options is required for subsonic inlet
    density=1.225;        # kg/m^3
    #pressure=101325;      # Pascal
    #temperature=300;      # Kelvin
}

Outlet

Outlet is a flexible type that will activate the proper boundary condition based on what is specified. Several different outlet implementations are available as detailed below.

BC_1 {
    type=outlet;

    # One of these options is required for subsonic outlet
    pressure=101325;      # Pascal
    #temperature=300;      # Kelvin
    #density=1.225;        # kg/m^3
    # If none of the above options are specified,
    # the outlet extrapolated everything from inside
    # This is only valid for supersonic outlets

    # Optional sponge layers to dampen reflections
    spongelayers=2;
}

Additional controls are offered through the kind options:

# Handle reverse flow by damping
BC_1 {type=outlet; kind=dampreverse}
# Prevent reverse flow
BC_1 {type=outlet; kind=noreverse}
# Force supersonic outlet conditions
BC_1 {type=outlet; kind=forcesupersonic}

Mach outlet

You can set either fixed local or average Mach number as outlet conditions:

BC_1 {
    type=outlet;
    mach=0.1;           # Target local Mach number
}

BC_1 {
    type=outlet;
    averagemach=0.1;    # Target average Mach number
}

Exact

Applies an exact solution from the specified exact flow option:

BC_1 {type=exact}

Nozzle

The nozzle boundary condition is mainly used for launch simulations where the boundary condition is applied to the nozzle exit. It also comes with a few physics options to improve stability and control the output. A basic boundary condition can be specified via the following:

    BC_1 {
        type                 = nozzle
        datafiletype         = 1
        components           = [1]
        filename             = nozzle_data.dat       
        ispecies = 0   
    }

where the nozzle boundary condition takes an input file filename which has several supported file types datafiletype described below. Like the above boundary conditions, they are applied to a component on the triangulation specified by components. In a multispecies simulation, you must also specify ispecies which corresponds to the 0 based index of the mixture in the fluid section for which species the nozzle will eject.

A complete list of all the options is provided below, along with default options where applicable:

Main Options:

Variable	Description
filename	(Required) The filename of the describing boundary condition file
datafiletype	(Required) The type of datafile. 0 is a steady (constant) ascii file. 1 is radially varying constant ascii file, 2 is unsteady but constant in space ascii file, 3 is unsteady radially varying binary, 4 is unsteady radially varying ascii file
ispecies	(Required for multispecies/multiphase) Specifies the component of mixture for the nozzle to eject.
nozzlecenter	(Optional) Specifies the coordinates of the nozzle center
nozzlenormal	(Optional) Specifies the normal for the boundary condition to be applied on
asciilinetimestart	(Required) Only used for datafiletype=4. Specifies the line number to start reading the file.
asciilinetimestop	(Required) Only used for datafiletype=4. Specifies the line number to stop reading the file.
asciilinetimeskip	(Optional) Only used for datafiletype=4. Specifies the stride at which to read lines from the file. Default is 1.
asciilinereadspecies	(Optional) Only used for datafiletype=4. Bool for whether this file contains species information as well. Default is false.
numberofradialpoints	(Required) Only used for datafiletype=4. Specifies the number of radial points in the file.

Note

The mixture state prescribed when asciilinereadspecies=true takes precedent over ispecies. In this case, ispecies is disregarded by the reader.

Ramp Options:

Variable	Description
rampsteps	(Optional) Will linearly warmup the nozzle based on rampsteps number of timesteps by taking a linear combination of the reference Q and the user specified conditions in filename. Default is 0 (no step ramping).
rampnton	(Optional) The timestep number to turn the ramping on. Nozzle will be off prior to this timestep. (Only used with rampsteps). Default is -1.
ramptimestart	(Optional) If rampsteps is not turned on, the other ramping option is using simulation time. This is the same linear warmup between ramptimestart and ramptimestop. Before ramptimestart the nozzle is off if it is specified. Default is no time ramping.
ramptimestop	(Optional) The simulation end time to ramping.
timeshift	(Optional) Timeshift can be used to shift the nozzles in time relative to the simulation time. This is compatible with ramping. Default is 0.

Note

One should only use (rampsteps and rampnton) or (ramptimestart and ramptimestop). They will not work in a combined fashion.*

Physics Options:

Variable	Description
multiphaseequilibrium	(Optional) Only used in cartesian + multiphase. This computes multiphase equilibrium on the boundary condition. Default is false.
dostagnation	(Optional) Apply stagnation inlet boundary conditions. Not supported for multiphase. Default is false.
doautomaticsonicinlet	(Optional) Only used when dostagnation = false. If the boundary condition goes below subsonic speed, we get pressure from the interior. This keeps the velocity at the inlet >= mach 1. Not supported for multiphase. Default is false.
dothroatsonicinlet	(Optional) Only used when dostagnation = true. Forces the pressure to be sonic + 1%. Not supported for multiphase. Default is false.

ASCII Data Format

datafiletype = 0 or 1: steady, radially constant or varying, ascii file: radius(m) Pressure(Pa) Temperature(K) u_normal(m/s) u_radial(m/s)

datafiletype = 2: unsteady, radially constant, ascii file: time(s) Pressure(Pa) Temperature(K) u_normal(m/s) u_radial(m/s)

Info

Lines starting with # are omitted by the parser for datafiletype = 0, 1, 2.

datafiletype = 4: unsteady, radially varying, ascii file: line_number time(s) radius_x(m) radius_y(m) [radius_z(m)] Pressure(Pa) u_normal(m/s) u_radial(m/s) [u_azimuthal(m/s)] Temperature (m) [Species Mass Fractions]

Note

The format of datafiletype = 4 depends on the dimension of the configuration. The extra fields radius_z(m) and u_azimuthal(m/s) are required when running with 3D and must be omitted when running with 2D. Currently, the field u_azimuthal(m/s) is disregarded in the nozzle bc, however, it must still be present in the nozzle data file. Also, the fiels Species Mass Fractions are only present when asciilinereadspecies = true. Also, note that no comments are allowed in this file type format.

An example of a type-4 nozzle data base is shown below for the configuration asciilinetimestart = 0; asciilinetimestop = 30; asciilinetimeskip = 10; asciilinereadspecies = true:

  0   -10000  0.000 0 0  10000  1000   0  0    1500  0. 0. 1. 0.
  0   -10000  0.500 0 0  40000  1000   0  0    1500  0. 0. 1. 0.
  0   -10000  1.000 0 0  80000  1000   0  0    1500  0. 0. 1. 0.
  0   -10000  1.965 0 0 110000  1000   0  0    1500  0. 0. 1. 0.
  10   1.000  0.000 0 0  10000  1000   0  0    1500  0. 0. 1. 0.
  10   1.000  0.500 0 0  40000  1000   0  0    1500  0. 0. 1. 0.
  10   1.000  1.000 0 0  80000  1000   0  0    1500  0. 0. 1. 0.
  10   1.000  1.965 0 0 110000  1000   0  0    1500  0. 0. 1. 0.
  20   10000  0.000 0 0  10000  1000   0  0    1500  0. 0. 1. 0.
  20   10000  0.500 0 0  40000  1000   0  0    1500  0. 0. 1. 0.
  20   10000  1.000 0 0  80000  1000   0  0    1500  0. 0. 1. 0.
  20   10000  1.965 0 0 110000  1000   0  0    1500  0. 0. 1. 0.
  30   20000  0.000 0 0  10000  1000   0  0    1500  0. 0. 1. 0.
  30   20000  0.500 0 0  40000  1000   0  0    1500  0. 0. 1. 0.
  30   20000  1.000 0 0  80000  1000   0  0    1500  0. 0. 1. 0.
  30   20000  1.965 0 0 110000  1000   0  0    1500  0. 0. 1. 0.

Distributed

Distributed boundary conditions use externally-defined flow data that is mapped to boundary faces. It uses surface solution files typically sampled from a previous solution. Currently cgns format is supported.

BC_1 {
    type=distributed;
    # Required: File containing distributed BC data
    filename=path/to/distributed_data.cgns;

    # Optional: Transform coordinates from data file to simulation
    scalecoords=[1.0,1.0,1.0];      # Scale coordinates [default=[1,1,1]]
    translatecoords=[0.0,0.0,0.0];  # Translate coordinates [default=[0,0,0]]
    rotatecoordsdeg=[0.0,0.0,0.0];  # Rotate coordinates in degrees [default=[0,0,0]]

    # Optional: Transform velocities from data file to simulation
    scalevelocity=[1.0,1.0,1.0];      # Scale velocity components [default=[1,1,1]]
    translatevelocity=[0.0,0.0,0.0];  # Translate velocities [default=[0,0,0]]
    rotatevelocitydeg=[0.0,0.0,0.0];  # Rotate velocities in degrees [default=[0,0,0]]
}

For unsteady distributed boundary conditions, use the following format:

BC_1 {
    type=distributed;
    isunsteady=true;

    # Required: Time step range and format specification
    ntlo=1000;       # First time step to use
    nthi=2000;       # Last time step to use
    ntskip=10;       # Interval between time steps
    numdigits=6;     # Number of digits in filename formatting
    prefix="data_";  # Prefix for filenames
    suffix=".cgns";    # Suffix for filenames

    # Optional: Transform specifications as in steady case
    scalecoords=[1.0,1.0,1.0];
    # ... other transform options
}

The distributed BC files should contain flow variables that are mapped to boundary faces using either nearest-point sampling or interpolation. For unsteady boundary conditions, the solver will interpolate between time steps as needed.

Note

Files for unsteady BCs are expected to follow the naming pattern: prefix + formatted_time_step + suffix For example: with prefix="data_", numdigits=6, suffix=".cgns", and time step 1000, the file would be "data_001000.cgns"

Warning

Distributed boundary conditions are computationally expensive and not supported when using GPU acceleration.

References

---

1. Sozer, E., Ghate, A. S., Kenway, G. K., Barad, M. F., Sousa, V. C., and Kiris, C. C., "Evaluation of Voronoi Meshes for Large Eddy Simulations of High Lift Aerodynamics," American Institute of Aeronautics; Astronautics, 2023. https://doi.org/10.2514/6.2023-0255 ↩

2. Sozer, E., Barad, M. F., Housman, J. A., Cadieux, F., Sousa, V., Dumlupinar, E., Saurabh, K., Ravikumar, K., and Duensing, J. C., "LAVA WMLES Results for the 5th High-Lift Prediction Workshop," American Institute of Aeronautics; Astronautics, 2025. https://doi.org/10.2514/6.2025-0887 ↩