[Eilmer Reference Manual, v4]

In a nutshell, transient flow solvers e4shared and e4mpi are gas-flow simulation codes. The set up of a simulation involves writing an input script that defines the domain of the gas flow as one or more FluidBlock objects with initial gas flow states and prescribed boundary conditions. Over small time steps, the code then updates the flow state in each cell within the flow domain according to the constraints of mass, momentum and energy, thus allowing the flow to evolve subject to the applied boundary conditions. The following sections provide brief details on many items that might go into your input script. Note that this document is for reference, after you have read the guides at https://gdtk.uqcloud.net/docs/eilmer/user-guide/ .

1. Configuration options

There are a large number of configuration options that can be set in the input script. The options are set in the input script by adding lines of the form:

config.option = value

Here are all of the available configuration options and the default values if left unset. Note you do not have to set all of these values in the input script. In a typical input script, one might set 10-15 input variables.

1.1. Geometry

config.dimensions

int, default: 2
Spatial dimensionality, 2 or 3

config.axisymmetric

bool, default: false
Sets axisymmetric mode for 2D geometries. The axis of rotation aligns with the x-axis (at y=0).

1.2. Time-stepping

config.gasdynamic_update_scheme

string, default: 'predictor-corrector'
The options are: 'euler', 'pc', 'predictor-corrector', 'midpoint', 'classic-rk3', 'tvd-rk3', 'denman-rk3', 'moving-grid-1-stage', 'moving-grid-2-stage', 'backward_euler', 'implicit_rk1'. Most of the update schemes are explicit, except for the last two listed here. Note that 'pc' is equivalent to 'predictor-corrector'. If you want time-accurate solutions, use a two- or three-stage explicit stepping scheme, otherwise, (explicit) Euler stepping has less computational expense but you may get less accuracy and the code will not be as robust for the same CFL value. For example the shock front in the Sod shock tube example is quite noisy for Euler stepping at CFL=0.85 but is quite neat with any of the two- or three-stage stepping schemes at the same value of CFL. The explicit midpoint and predictor-corrector schemes produce a tidy shock up to CFL = 1.0 and the rk3 schemes still look tidy up to CFL = 1.2. If you care mainly about the end result, you might be able to use an implicit scheme with larger CFL number. We have run the 'backward_euler' scheme up to a CFL value of 50 for some of the turbulent flat-plate examples. The implicit schemes are more expensive than the low-order explicit schemes, per time-step and, at large CFL numbers, will not be time accurate. The trade-off will be something that you have to explore but, if you care only about a steady state, the 'backward_euler' scheme may get you there more quickly than any of the explicit schemes. Your mileage will vary. Note that you may use dashes or underscores when spelling the scheme names.

config.fixed_time_step

boolean, default: false
Normally, we allow the time step size to be determined from cell conditions and CFL number. Setting config.fixed_time_step=true forces the time step size to be unchanged from config.dt_init.

config.cfl_value

float, default: 0.5
The CFL number is ratio of the time-step divided by the time for the fastest signal to cross a cell. The time-step adjustment process tries to set the time-step for the overall simulation to a value such that the maximum CFL number for any cell is at this cfl_value. If you are having trouble with a simulation that has lots of sudden flow field changes, decreasing the size of cfl_value may help.

config.cfl_schedule

table of float pairs, default {{0.0, config.cfl_value},}
If you want to vary the cfl_value as the simulation proceeds, possibly starting with small values and then increasing them as the flow field settles to something not too difficult for the flow solver, you can specify a table of {time,cfl_value} pairs that are interpolated by the solver. Outside the specified range of times, the cfl_value will be the closest end-point value.

config.dt_init

float, default: 1.0e-3
Although the computation of the time step size is usually automatic, there might be cases where this process does not select a small enough value to get the simulation started stably. For the initial step, the user may override the computed value of time step by assigning a suitably small value to config.dt_init. This will then be the initial time step (in seconds) that will be used for the first few steps of the simulation process. Be careful to set a value small enough for the time-stepping to be stable. Since the time stepping is synchronous across all parts of the flow domain, this time step size should be smaller than half of the smallest time for a signal (pressure wave) to cross any cell in the flow domain. If you are sure that your geometric and boundary descriptions are correct but your simulation fails for no clear reason, try setting the initial time step to a very small value. For some simulations of viscous hypersonic flow on fine grids, it is not unusual to require time steps to be as small as a nanosecond.

config.dt_max

float, default: 1.0e-3
Maximum allowable time step (in seconds). Sometimes, especially when strong source terms are at play, the CFL-based time-step determination does not suitably limit the size of the allowable time step. This parameter allows the user to limit the maximum time step directly.

config.viscous_signal_factor

float, default: 1.0
By default, the full viscous effect for the signal calculation will be used within the time-step calculation. It has been suggested that the full viscous effect may not be needed to ensure stable calculations for highly-resolved viscous calculations. A value of 0.0 will completely suppress the viscous contribution to the signal speed calculation but you may end up with unstable stepping. It’s a matter of try a value and see if you get a larger time-step while retaining a stable simulation.

config.stringent_cfl

boolean, default: false
The default action for structured grids is to use different cell widths in each index direction. Setting config.stringent_cfl=true will use the smallest cross-cell distance in the time-step size check.

config.cfl_count

int, default: 10
The number of time steps between checks of the size of the time step. This check is expensive so we don’t want to do it too frequently but, then, we have to be careful that the flow does not develop suddenly and the time step become unstable.

config.max_time

float, default: 1.0e-3
The simulation will be terminated on reaching this value of time.

config.max_step

int, default: 100
The simulation will be terminated on reaching this number of time steps. You will almost certainly want to use a larger value, however, a small value is a good way to test the starting process of a simulation, to see that all seems to be in order.

config.dt_plot

float, default: 1.0e-3
The whole flow solution will be written to disk when this amount of simulation time has elapsed, and then again, each occasion the same increment of simulation time has elapsed.

config.dt_history

float, default: 1.0e-3
The history-point data will be written to disk repeatedly, each time this increment of simulation time has elapsed. To obtain history data, you will also need to specify one or more history points.

config.with_local_time_stepping

boolean, default: false
If a steady-state solution is sought, then it is possible to accelerate the explicit time-stepping scheme by allowing each cell in the domain to integrate in time with their own time-step that satisfies only the local CFL condition. This is sometimes referred to as Local Time-Stepping (LTS). Setting config.with_local_time_stepping=true enables this mode of time-stepping. Note that this mode of time-stepping is only suitable for steady-state simulations.

config.local_time_stepping_limit_factor

int, default: 10000
To improve the stability of simulations that enable config.with_local_time_stepping, users can limit the largest allowable local time-step as a function of the smallest global time-step. The allowable time-step for each cell is set to be the minimum of the cell’s local time-step and the global minimum time-step multiplied by config.local_time_stepping_limit_factor. The config.local_time_stepping_limit_factor is typically set to a value between 500 and 10000.

config.residual_smoothing

boolean, default: false
Residual smoothing essentially applies an averaging filter on the residuals during the time-integration procedure. Residual smoothing has two advantages: (1) The averaging increases the allowable CFL for a particular explicit time-stepping scheme; and (2) The averaging ensures that the local time-step for two adjacent cells doesn’t differ by a large magnitude when using config.with_local_time_stepping, this has a stabilising effect.

config.residual_smoothing_type

string, default: 'explicit'
The options are: 'explicit', 'implicit'. When config.residual_smoothing=true, the residuals can either be smoothed via an explicit averaging method or an implicit averaging method. Explicit averaging has less computational expense, however, implicit averaging typically results in a larger stable time step.

config.residual_smoothing_weight

float, default: 0.2
A weighting factor used in the residual averaging when config.residual_smoothing=true.

config.residual_smoothing_iterations

int, default: 2
The number of Jacobi iterations used for the implicit residual averaging when config.residual_smoothing=true and config.residual_smoothing_type='implicit'.

1.3. Block marching

config.block_marching

boolean, default: false
Normal time iteration proceeds on all blocks simultaneously, however, such a time-marching calulation may be very expensive computationally. Setting config.block_marching=true enables a sequencing of the time integration such that at any one instant, only two slices of blocks are being integrated. The i-direction is the marching direction and the assumed dominant (supersonic) flow direction. The blocks are assumed to be in a regular array with a fixed number of blocks in the j- and k-directions to the entire flow domain.

config.nib, config.njb, config.nkb

int default: 1, 1, 1
are the number of blocks in each index direction. To make the best of block marching, you should have config.nib set to a fairly large number. Since the array of blocks is assumed regular, you cannot have very complicated geometries. Simple ducts, nozzles and plates are the intended applications. As seen in the examples, it may be convenient to define the full domain with one or more calls to the constructor FBArray:new. There is a restriction that the overall flow domain be assembled as a single structured array of FlowBlock objects.

config.propagate_inflow_data

boolean, default: false
By default, the integration begins in each set of blocks from the initial gas state set up in the preparation phase of the simulation. Some advantage may be gained following integration of the first block slices by initializing subsequent block slices with the downstream (east boundary) flow states. Setting config.propagate_inflow_data=true propagates these data across each new block slice, before the integration process for the slice begins.

config.save_intermediate_results

boolean, default: false
Usually, a single set of solution files (after marching over all block slices) is all that is required. Sometimes, when debugging a troublesome calculation, it may be useful to have a solution written after the time-integration process for each pair of block slices. Set this parameter true to get these intermediate solutions written.

1.4. Spatial reconstruction

config.interpolation_order

int, default: 2
Before applying the flux calculator, high-order reconstruction is applied. Setting config.interpolation_order=1 results in no reconstruction of intra-cell flow properties.

config.interpolation_delay

float, default: 0.0
The time (in seconds) to wait before increasing the interpolation order. If the delay is enabled, the simulation will begin using interpolation_order=1. After the delay time, the interpolation order will switch to that specified in the configuration. This delay is useful to help with robustness during transient starts because the first order spatial reconstruction is very reliable (albeit less accurate as compared to high order reconstruction).

config.apply_limiter

boolean, default: true
By default, we apply a limiter to the flow-field reconstruction.

config.extrema_clipping

boolean, default: true
By default, we do extrema clipping at end of each scalar-field reconstruction. Setting config.extrema_clipping=false suppresses clipping.

config.thermo_interpolator

string, default: 'rhou'
String to choose the set of interpolation variables to use in the interpolation, options are 'rhou', 'rhop', 'rhoT' and 'pT'.

1.5. Flux calculator

config.flux_calculator

string, default: 'adaptive_hanel_ausmdv'
Selects the flavour of the flux calculator. Options are:

  • 'efm' A cheap and very diffusive scheme by Pullin and Macrossan. For most hypersonic flows, it is too diffusive to be used for the whole flow field but it does work very nicely in conjunction with AUSMDV, especially for example, in the shock layer of a blunt-body flow.

  • 'ausmdv' A good all-round scheme with low-diffusion for supersonic flows.

  • 'adaptive_efm_ausmdv' A blend of the low-dissipation AUSMDV scheme for the regions away from shocks with the much more diffusive EFM used for cell interfaces near shocks. It seems to work quite reliably for hypersonic flows that are a mix of very strong shocks with mixed regions of subsonic and supersonic flow. The blend is controlled by the parameters config.compression_tolerance and config.shear_tolerance that are described below.

  • 'ausm_plus_up' Implemented from the description by MS Liou (2006). It should be accurate and robust for all speed regimes. It is the flux calculator of choice for very low Mach number flows, where the fluid behaviour approaches the incompressible limit. For best results, you should set the value of M_inf.

  • 'hlle' The Harten-Lax-vanLeer-Einfeldt (HLLE) scheme. It is somewhat dissipative and is the only scheme usable with MHD terms.

  • 'adaptive_hlle_ausmdv' As for 'adaptive_efm_ausmdv' but with the dissipative scheme being the HLLE flux calculator.

  • 'hanel' The Hanel-Schwane-Seider scheme, from their 1987 paper. It also dissipative and is somewhat better behaved than our EFM implementation.

  • 'adaptive_hanel_ausmdv' As for 'adaptive_efm_ausmdv' but with the dissipative scheme being the Hanel-Schwane-Seider flux calculator.

  • 'roe' The Phil Roe’s classic linearized flux calculator.

  • 'adaptive_hlle_roe' A blend of Roe’s low-dissipation scheme and the more dissipative HLLE flux calculator.

  • 'adaptive_ausmdv_asf' A blend of Summation-by-Parts method of Fisher (high order, low-dissipation scheme) and the ausmdv method.

The default adaptive scheme is a good all-round scheme that uses AUSMDV away from shocks and Hanel-Schwane-Seider flux calculator near shocks.

Any of the schemes with the prefix 'adaptive' can be used in a blended mode if the shock_detector_smoothing is enabled.

config.compression_tolerance

float, default: -0.30
The value of relative velocity change (normalised by local sound-speed) across a cell-interface that triggers the shock-point detector. A negative value indicates a compression. When an adaptive flux calculator is used and the shock detector is triggered, the more-dissipative flux calculation will be used in place of the default low-dissipation calculation. A value of -0.05 seems OK for the Sod shock tube and sharp-cone inviscid flow simulations, however, a higher value is needed for cases with viscous boundary layers, where it is important to not have too much numerical dissipation in the boundary layer region.

config.shear_tolerance

float, default: 0.20
The value of the relative tangential-velocity change (normalised by local sound speed) across a cell-interface that suppresses the use of the high-dissipation flux calculator even if the shock detector indicates that high-dissipation scheme should be used within the adaptive flux calculator. The default value is experimentally set at 0.20 to get smooth shocks in the stagnation region of bluff bodies. A smaller value (say, 0.05) may be needed to get strongly expanding flows to behave when regions of shear are also present.

config.M_inf

float, default: 0.01
representative Mach number for the free stream. Used by the ausm_plus_up flux calculator.

config.shock_detector_smoothing

int, default = 0
How many iterations of shock detector averaging. Higher values spread the influence of the shock detector further from the shock, increasing numerical dissipation but providing a smoother transition between flux calculators in adaptive schemes.

Tip
You can read more about this in the technical note Hybrid fluxes, shock detection and smoothing.
config.strict_shock_detector

boolean, default: true
If this is true then any face that is part of a cell with S>0 will have its own S value set to 1. This is to emulate legacy behaviour and thus defaults to true. For true blending of fluxes this should be set to false.

config.frozen_shock_detector

boolean, default: false
If set to true, the shock detector field will be frozen after the number of time steps prescribed by config.shock_detector_freeze_step. Used to prevent flickering of adaptive flux methods between their constituent flux calculators.

config.shock_detector_freeze_step

int, default: 1000
The time step at which the shock detector field will be frozen.

1.6. Viscous effects

config.viscous

boolean, default: false
If set true, viscous effects will be included in the simulation.

config.separate_update_for_viscous_terms

boolean, default: false
If set true, the update for the viscous transport terms is done separately to the update for the convective terms. By default the updates are done together in the gas-dynamic update procedure.

config.viscous_delay

float, default: 0.0
The time (in seconds) to wait before applying the viscous terms. This might come in handy when trying to start blunt-body simulations.

config.viscous_factor_increment

float, default: 0.01
The per-time-step increment of the viscous effects, once simulation time exceeds config.viscous_delay.

config.turbulence_model

string, default: 'none'
String specifying which model to use. Options are: 'none' 'k_omega' 'spalart_allmaras'

config.turbulence_prandtl_number

float, default: 0.89

config.turbulence_schmidt_number

float, default: 0.75

config.max_mu_t_factor

float, default: 300
The turbulent viscosity is limited to laminar viscosity multiplied by this factor.

config.transient_mu_t_factor

float, default: 1.0

1.7. Finite-rate thermo-chemistry

config.reacting

boolean, default: false
Set to true to activate the finite-rate chemical reactions.

config.reactions_file

string, default: 'chemistry.lua'
File name for reaction scheme configuration.

config.reaction_time_delay

float, default: 0.0
Time after which finite-rate reactions are allowed to start.

config.reaction_fraction_schedule

table of float pairs, default {{0.0, 1.0},}
If you want to vary the fraction of the reaction time step per gas-dynamic time step as the simulation proceeds, you can specify a table of {time, fraction} pairs that are interpolated by the solver. If the initial flow is likely to have a harsh shock reflection with high temperatures and correspondingly high reaction rates, it may be beneficial to discard time accuracy and start with a small fraction, increasing it (to a maximum of 1.0) as the flow field settles. Outside the specified range of times, the fraction will be the closest end-point value.

config.T_frozen

float, default: 300.0
Temperature (in degrees K) below which reactions are frozen. The default value is 300.0 since most reaction schemes seem to be valid for temperatures above this, however, you may have good reasons to set it higher or lower.

1.8. Species Diffusion

config.mass_diffusion_model

string, default: 'none'
Set the expression for laminar diffusion flux of species along concentration gradients. Valid options are: 'none' (no diffusion) and 'ficks_first_law' (self explanatory). Note that this option is ignored in turbulent flow, where the turbulent diffusivity is used instead.

config.diffusion_coefficient_type

string, default: 'none'
Set the type of formula used to compute the species diffusion coefficients in a diffusion enabled simulation. Valid options are: 'constant_lewis_number' (works with any multi-species model), 'species_specific_lewis_numbers' (requires specific gas models), and 'binary_diffusion' (requires specific gas models). Note that leaving this parameter as 'none' in a simulation where config.mass_diffusion_model='ficks_first_law' will cause a run-time error.

config.lewis_number

float, default: 1.0
Set the value for the Lewis number, which computes the species diffusion coefficients based on the thermal conductivity, in a diffusion enabled simulation with config.diffusion_coefficient_type=constant_lewis_number.

1.9. Special initialisation

config.diffuse_wall_bcs_on_init

bool, default: false
Set to diffuse(/blend) conditions at the wall into the flow domain as an initial condition

config.number_init_passes

int, default: 30
Set how many passes the diffusion is applied. Each pass diffuses another layer of cells into domain. So 10 passes would effect 10 layers of cells from the wall (in a structured grid).

config.wall_temperature_on_init

float, default: -1.0
Set the wall temperature to use when diffusing the wall conditions. This value must be set when an adiabatic wall condition is selected since there is no initial guess of the wall temperature. If this value is set for a fixed temperature wall, then this value overrides the wall temperature for the purposes of this diffusion-style initialisation. It does not change the wall temperature selected in the boundary condition. A value of -1.0 indicates that this selection is not active and the wall boundary temperature should be used. This is default action.

The following lines are an example of using the diffusion-style initialisation. We assume an adiabatic wall and give a starting guess of 600.0 K for the wall temperature. The blending operation is performed 20 times in this example:

config.diffuse_wall_bcs_on_init = true
config.number_init_passes = 20
config.wall_temperature_on_init = 600.0

1.10. Miscellaneous

config.title

string, default: "Eilmer4 simulation"
Title for the simulation

config.adjust_invalid_cell_data

boolean, default: false
Usually, you will want the flow solver to provide its best estimate for your flow, however, there are flow situations for which the flow solver will not compute physically valid flow data. If you encounter a difficult flow situation and are prepared to fudge over a few cells, then set this parameter to true and max_invalid_cells to a non-zero value. Be cautious when using this option and use it only when you have exhausted more reasoned options. If there is a problem that is more than just a difficult patch of flow that will blow by, it may allow you to go further into a bad situation and get even more confused about what the underlying issue really is.

config.max_invalid_cells

int, default: 0
The maximum number of bad cells that will be tolerated on decoding conserved quantities. If this number is exceeded, the simulation will stop.

config.report_invalid_cells

boolean, default: true
If you are stuck with having to fudge over cells, you probably will want to know about them until, of course, that you don’t. Set this parameter to false to silence the reports of bad cells being fudged over.

config.apply_bcs_in_parallel

boolean, default, true
This will be the fastest calculation, however, some boundary conditions, such as the shock-fitting need to cooperate across blocks and so will have race conditions if applied in parallel. If your simulation has such a boundary condition, set this parameter to false to favour safety above speed.

config.udf_source_terms

boolean, default: false
Set to true to apply user-defined source terms, as supplied in a Lua file.

config.udf_source_terms_file

string, default: "dummy-source-terms.txt"
Name of the Lua file for the user-defined source terms.

config.print_count

int, default: 20
Number of time steps between printing status information to the console.

config.control_count

int, default: 10
Number of time steps between re-parsing the job`.control` file. If the job`.control` has been edited, then the new values are used after re-parsing.

config.MHD

boolean, default: false
Set to true to make MHD physics active.

config.do_temporal_DFT

boolean, default: false
Set to true to calculate the temporal Discrete Fourier Transform of the pressure on a cell-by-cell basis on-the-go. Useful for large jobs where storing the number of snapshots required for Fourier decomposition is prohibitive. Results written to file DFT/jobName.DFT.bxxxx.tar.gz, with a line containing the real and imaginary components of each Fourier mode i.e. Re(A0), Im(A0), Re(A1), Im(A1), …​ Im(A_n). Requires setting of additional config options config.DFT_n_modes and config.DFT_step_interval. Use of config.fixed_time_step is also very strongly recommended, non-fixed time stepping will lead to erronous mode amplitudes and phases.

config.DFT_n_modes

int, default: 5
The number of Fourier modes to be computed in the temporal Discrete Fourier Transform. Should be equal to config.max_steps / config.DFT_step_interval (see next entry).

config.DFT_step_interval

int, default: 10
The number of steps between each increment of the temporal Discrete Fourier Transform modes. Should be equal to config.max_steps / config.DFT_n_modes.

1.11. Special Zones

Sometimes you may wish to apply specific configuration only for a limited region. There are special zones that can be constructed for turbulence, reactions and ignition and order of reconstruction. The region of application is defined by a pair of diagonally-opposing points that specify the minimum and maximum coordinate values.

TurbulentZone:new{p0={x=x0, y=y0}, p1={x=x1, y=y1}}

The turbulence model is active throughout the flow but its effect on the flow field is masked outside of any defined turbulent zones. This is achieved by the code setting the turbulence viscosity and conductivity to zero for finite-volume cells that fall outside of all regions defined as a TurbulentZone. If there are no such defined regions, the whole flow field is allowed to have nonzero turbulence viscosity.

ReactionZone:new{p0={x=x0, y=y0}, p1={x=x1, y=y1}}

In a flow with an active reaction scheme, this type of zone makes it possible to selectively allow reactions to proceed, or not. If the centre of a cell lies within the reaction zone, the finite-rate chemistry is allowed to proceed, else the species fractions are maintained constant. If no reaction zones are specified and a reaction scheme is active, then reactions are permitted for the entire flow field.

IgnitionZone:new{p0={x=x0, y=y0}, p1={x=x1, y=y1}, T=Tig}

Tig controls the reaction rate used for chemical reactions, without effecting the gas temperature in the flow field. The rate-controlling temperature is used to evaluate the chemical reaction rates only within the physical extents of the ignition zone. The effect of this zone can be limited in time by specifying a nonzero values for config.ignition_time_start and config.ignition_time_stop. While the zone is active, the reaction rates within the zone are altered. The rate-controlling temperature is typically set to an artificially inflated value to promote ignition.

SuppressReconstructionZone:new{p0={x=x0, y=y0}, p1={x=x1, y=y1}}

Inside such a zone, high-order reconstruction of the flow properties is suppressed. This can sometimes help with difficult flow features, such as an expansion around a sharp corner.

2. FlowState objects

Constructing a FlowState object is a way to collect the full set of data that is required for defining the flow state initially throughout a block and also for boundaries that have inflow.

FlowState:new{p, T, T_modes, massf, velx, vely, velz, tke, omega, mu_t, k_t}

FlowState objects are always constructed in the context of a gas model that defines the set of chemical species and other parameters of the gas. For a quiescent, single-species gas the minimum data that you need to specify are:

p

float, no default
Gas pressure, in Pa.

T

float, no default
Gas temperature, in K.

If the flow has nonzero velocity, you can specify the velocity components:

velx

float, default 0.0
Velocity in x-direction, m/s.

vely

float, default 0.0
Velocity in y-direction, m/s.

velz

float, default 0.0
Velocity in z-direction, m/s.

If the gas is composed of more than one chemical species, you will need to specify the mass-fractions of those species.

massf

table of named float values
Mass fractions of the species specified as a table with named mass-fraction values. The default is a table with the first species at mass-fraction 1.0 and all others at 0.0. For example, the ideal-air gas model has a default {'air'=1.0,}. You need specify only the species with nonzero mass fractions, but that collection of specified values must add to 1.0.

T_modes

table of float values
If the gas model has modes of internal energy storage beyond the usual thermal modes, you may specify the temperatures corresponding to those energy modes. If you do not provide a table, the default value for all modes will be T.

When a turbulence model has been selected, you may also specify the local intensity as:

tke

float
Turbulent kinetic energy.

omega

float
Turbulent frequency.

mu_t

float
Turbulent viscosity.

k_t

float
Turbulent conductivity.

Once a FlowState object has been constructed, you have access to the other GasState data via the usual Lua table access. For example you could access the local speed of sound and compute a mach number

inflow = FlowState:new{p=5.0e3, T=300.0, velx=1000.0}
mach = inflow.velx / inflow.a

3. Block objects

3.1. FluidBlock objects

We decompose the fluid simulation domain into blocks. Our simplest decomposition would be one block to represent the entire domain). We call these units of decomposition FluidBlocks. The input script is largely focussed on configuring FluidBlocks to represent the fluid simulation domain. The configuration includes setting boundary conditions and the initial state. The full list of options are described here.

FluidBlock:new{grid, initialState,
               bcList|bcDict,
               active, label, omegaz, may_be_turbulent}
grid

StructuredGrid or UnstructuredGrid object, no default
A grid object must be supplied. The type of grid (structured or unstructured) sets the type of block. The grid dimensions should also be compatible with the global setting for dimenions. (See: config.dimensions)

initialState

FlowState, FlowSolution or a user-defined function
The fluid simulation domain requires an initial set of flow properties throughout the domain. These flow properties are set block-by-block using this initialState parameter. There is no default: the user must supply something for this parameter.

Most commonly, a FlowState object is given as the initial state for a block. This sets spatially-uniform flow properties across the entire block.

A FlowSolution object might also be used to set the initial flow state in a block. Usually the FlowSolution comes from a previously completed simulation. One use for this is to start a finer resolution calculation using a coarser simulation. Another use could be to do a staged simulation: when extending the domain in a staged simulation, one might want to start with some blocks initialised from an earlier truncated domain simulation. The cells in the block are initialised by searching the FlowSolution for the nearest corresponding cell in the old solution and copying that value. No attempt at interpolating based on nearby cells or mismatch in cell centres is made.

A user-defined function can be provided as a Lua function. This allows one to provide a spatially-varying initial state. The function contract is: accepts x, y and z as parameters; returns a FlowState object. An example of initialising a shock tube with high and low pressure sections is shown in the 3D Sod shock tube case. See: dgd/examples/eilmer/3D/sod-shock-tube/sg/sod.lua. That initialisation function is repeated here:

function tube_gas(x, y, z)
   if x < 0.5 then
      p = 1.0e5; T = 348.4
   else
      p = 1.0e4; T = 278.8
   end
return FlowState:new{p=p, velx=0.0, vely=0.0, T=T}
end
bcList

table: {north=BCObject, east=BCObject, south=BCObject, west=BCObject, [top=BCObject, bottom=BCObject]}
The bcList is used to set boundary conditions on the faces of blocks with an underlying structured grid. Each table key accepts a boundary condition object. If any table entry is missing, the default-supplied boundary condition is a slip wall, WallBC_WithSlip. In 3D, the top and bottom faces can also be set. If a boundary is a connection, this can be explicitly set or you can rely on the identifyBlockConnections() helper function to make that connection for you. That helper function is called after the blocks have been configured. An example to set a supersonic inflow on a west boundary, and an outflow on the north boundary is:

bcList={west=InFlowBC_SupersonicIn:new{flowState=inflow},
        north=OutFlowBC_Simple:new{}}
bcDict

table: {tagA=BCObject, tagB=BCObject, …​}
The bcDict is used to set boundary conditions on faces of blocks with an underlying unstructured grid. The tag strings in the table should match the tags used in the underlying unstructured grid. If a boundary face is not set explicitly in this table, then a default boundary condition of a slip wall is applied (WallBC_WithSlip).

active

bool, default: true
By default, blocks are marked as active for inclusion in the simulation. You could disable a block’s participation in the simulation domain, setting this parameter to false. Be careful of how this block is connected to others if setting it to inactive; the flow state will not be updated in an inactive block.

label

string, default: FluidBlock-<id>
Set this parameter to give the block a label. This is a convenience for the user and might be useful for some custom post-processing to work with blocks that have a specific label.

omegaz

float, default: 0.0
angular speed (in rad/s about the z-axis) for the rotating frame of reference.
Useful for the calculation of flow in turbomachines. The default value of 0.0 implies no rotating frame.

may_be_turbulent

bool, default: true
Normally, when a turbulence model is active, turbulent flow will be allowed in all blocks. This parameter may be set false to suppress the turbulence actions for the block. Use sparingly.

3.2. FluidBlock arrays

To generate more than one FluidBlock on a single, structured grid, use

FBArray:new{grid, nib, njb, nkb,
            initialState,
            bcList|bcDict,
            active, label, omegaz, may_be_turbulent}

The entries to the table are mostly the same as for the construction of a single FluidBlock object. The extra arguments are:

nib

int, default 1
The number of blocks in the i-index direction.

njb

int, default 1
The number of blocks in the j-index direction.

nkb

int, default 1
The number of blocks in the k-index direction.

The table that is returned contains an arrays of all of the FluidBlock objects defined over the partitioned grid:

blockArray

will be a multi-dimensional array, indexed as [ib][jb][kb], with 1⇐ib⇐nib, 1⇐jb⇐njb, 1⇐kb⇐nkb.

blockCollection

will be a single-dimensional array, also starting at 1. Note that the index of the blocks within this array follow the Lua convention and will be different to the id assigned to the blocks, which will start at 0, and thus align with the D-language indexing.

3.3. SolidBlock objects

We decompose the solid simulation domain into blocks. We call these units of decomposition SolidBlocks. The configuration includes setting the solid properties, boundary conditions and the initial temperature. The full list of options are described here.

SolidBlock:new{grid, initTemperature, bcList, active, label, properties}

grid

StructuredGrid object, no default
A structured grid object must be supplied. The grid dimensions should also be compatible with the global setting for dimenions. (See: config.dimensions)

initTemperature

float, no default
The solid simulation domain requires an initial temperature throughout the domain. The temperature is set block-by-block using this initTemperature parameter. There is no default: the user must supply something for this parameter.

bcList

table: {north=BCObject, east=BCObject, south=BCObject, west=BCObject, [top=BCObject, bottom=BCObject]}
The bcList is used to set boundary conditions on the faces of blocks with an underlying structured grid. Each table key accepts a boundary condition object. If any table entry is missing, the default-supplied boundary condition is an adiabatic wall, SolidAdiabaticBC. In 3D, the top and bottom faces can also be set. If a boundary is a connection, this can be explicitly set or you can rely on the identifyBlockConnections() helper function to make that connection for you. That helper function is called after the blocks have been configured. An example to set a fixed temperature on a west boundary, and a constant flux on the north boundary is:

bcList={west=SolidFixedTBC:new{Twall=300.0},
        north=SolidConstantFluxBC:new{}}
active

bool, default: true
By default, blocks are marked as active for inclusion in the simulation. You could disable a block’s participation in the simulation domain, setting this parameter to false. Be careful of how this block is connected to others if setting it to inactive; the flow state will not be updated in an inactive block.

label

string, default: SolidBlock-<id>
Set this parameter to give the block a label. This is a convenience for the user and might be useful for some custom post-processing to work with blocks that have a specific label.

properties

table: {rho=rho, k=k, Cp=Cp}
The properties table is used to set the material-specific solid proprties where rho is the density, k is the thermal conductivity, and Cp is the specific heat capacity.

3.4. SolidBlock arrays

To generate more than one SolidBlock on a single, structured grid, use

SolidBlockArray:new{grid, nib, njb, nkb,
                    initTemperature,
                    bcList,
                    active, label, properties}

The entries to the table are mostly the same as for the construction of a single SolidBlock object. The extra arguments are:

nib

int, default 1
The number of blocks in the i-index direction.

njb

int, default 1
The number of blocks in the j-index direction.

nkb

int, default 1
The number of blocks in the k-index direction.

The table that is returned contains an arrays of all of the SolidBlock objects defined over the partitioned grid:

blockArray

will be a multi-dimensional array, indexed as [ib][jb][kb], with 1⇐ib⇐nib, 1⇐jb⇐njb, 1⇐kb⇐nkb.

blockCollection

will be a single-dimensional array, also starting at 1. Note that the index of the blocks within this array follow the Lua convention and will be different to the id assigned to the blocks, which will start at 0, and thus align with the D-language indexing.

4. Boundary Conditions

To complete a definition of a FluidBlock, boundary conditions for all block boundaries, need to be specified. These may be given as a list to the FluidBock constructor or they may be attached to particular boundaries of an already existing FluidBlock object. In the absence of your specification, the default boundary condition is a slip wall, WallBC_WithSlip.

4.1. Walls

4.1.1. Slip wall

WallBC_WithSlip:new{label, group}
or
WallBC_WithSlip0:new{label, group}

is used where we want a solid wall with no viscous effects. The effect of the solid wall is achieved with the use of ghost cells and an appropriate reflection of the velocity vector. This is the default boundary condition where no other condition is specified.

label

string, default ""
A an optional tag string that may be applied to the boundary.

group

string, default ""
A an optional tag string that identifies the group to which this boundary belongs. It can be used, at run time and in postprocessing, to group the data for several boundaries together. This might be for force and moment estimation, for example.

There are variants to the basic wall condition.

WallBC_WithSlip1:new{label, group}

achieves a similar effect to WallBC_WithSlip0 but without the use of ghost cell data. To do this, the code makes use of one-sided flux-calculators at the boundary. You will need to specify this variant of the boundary condition for a moving-grid simulation where this boundary wall may have non-zero velocity. The one-sided flux calculators will correctly handle the normal component of the wall velocity while, effectively, the tangential velocity will be ignored.

WallBC_WithSlip2:new{label, group}

is Ingo’s variant of WallBC_WithSlip for moving mesh simulations for walls with normal velocity. It uses ghost cells with appropriate velocity settings.

4.1.2. No-slip, fixed-temperature

WallBC_NoSlip_FixedT:new{Twall, wall_function, catalytic_type, wall_massf_composition, label, group}
or
WallBC_NoSlip_FixedT0:new{Twall, wall_function, catalytic_type, wall_massf_composition, label, group}

is used where we want convective and viscous effects to impose a no-slip velocity condition and a fixed wall temperature. The convective effect of the solid wall is achieved with the use of ghost cells and an appropriate reflection of the velocity vector. We need to set config.viscous=true to make the temperature-setting part of this boundary condition effective.

Twall

The fixed wall temperature, in degrees K.

wall_function

an optional analytical model for the turbulent wall stress and heat transfer, allowing the user to scrimp on small grid cells close to the wall. Suitable for use with the k-omega turbulence model only.

catalytic_type

an optional boundary condition used for species transport equations in a reacting simulation. Options are: none which assumes no catalytic interaction with the wall, equilibrium which assumes instantly completed reactions at the wall, and fixed_composition which sets the mass fractions to a user-defined value.

catalytic_type

A table of species mass fractions required for catalytic_type=fixed_composition.

`WallBC_NoSlip_FixedT1`

is variant of the boundary condition that does not use ghost cells and is suitable for moving-grid simulations.

4.1.3. No-slip, user-defined temperature

WallBC_NoSlip_UserDefinedT:new{Twall, wall_function, catalytic_type, wall_massf_composition, label, group}

is similar to the fixed-temperature wall but allows an arbitrary temperature function to be specified via a lua script. This may be useful for simulating a heated experimental model, where some specific part of the model is heated by electrical resistance heating to a known elevatated temperature. Optional arguments are the same as the fixed temperature case.

Twall

The name of the lua file used to specify the wall temperature. Inside the code, this file is passed to a UserDefinedInterface effect, which takes care of calling it once per face in the boundary, for each timestep. For further information, see the section regarding user defined interface effects in the PDF User Guide.

4.1.4. No-slip, adiabatic

WallBC_NoSlip_Adiabatic:new{wall_function, catalytic_type, wall_massf_composition, label, group}
or
WallBC_NoSlip_Adiabatic0:new{wall_function, catalytic_type, wall_massf_composition, label, group}

is used where we want convective and viscous effects to impose no-slip at the wall but where there is no heat transfer. Optional arguments are the same as the fixed temperature case. We need to set config.viscous=true to make this boundary condition effective.

WallBC_NoSlip_Adiabatic1:new{label, group}

is variant of the boundary condition that does not use ghost cells and is suitable for moving-grid simulations.

4.1.5. Translating-surface, fixed-temperature

WallBC_TranslatingSurface_FixedT:new{Twall, v_trans, label, group}

is used where we want convective and viscous effects to impose a specified translating velocity condition and a fixed wall temperature, in a simulation with no grid motion. By translating we mean that the (flat) wall is moving tangential to the block boundary such that the gas velocity is nonzero and aligned with the edge of the (fixed) grid. We need to set config.viscous=true to make this boundary condition fully effective. An example of use is Couette flow between moving plates.

Twall

The fixed wall temperature, in degrees K.

v_trans

Vector3, default: {x=0.0, y=0.0, z=0.0}
Vector velocity of the translating wall. The value may be specified as a table of three named (x,y,z) components.

4.1.6. Translating-surface, adiabatic

WallBC_TranslatingSurface_Adiabatic:new{v_trans, label, group}

is used where we want convective and viscous effects to impose a specified translating velocity condition but no heat transfer, in a simulation with no grid motion. By translating we mean that the (flat) wall is moving tangential to the block boundary such that the gas velocity is nonzero and aligned with the edge of the (fixed) grid. We need to set config.viscous=true to make this boundary condition fully effective.

v_trans

Vector3, default: {x=0.0, y=0.0, z=0.0}
Vector velocity of the translating wall. The value may be specified as a table of three named (x,y,z) components.

4.1.7. Rotating-surface, fixed-temperature

WallBC_RotatingSurface_FixedT:new{Twall, r_omega, centre, label, group}

is used where we want convective and viscous effects to impose a specified velocity condition and a fixed wall temperature on a circular or cylindrical fixed surface. By rotating we mean that the (curved) wall is moving tangential to the block boundary such that the gas velocity is nonzero and aligned the edge of the (fixed) grid. We need to set config.viscous=true to make this boundary condition fully effective.

Twall

The fixed wall temperature, in degrees K.

r_omega

Vector3, default: {x=0.0, y=0.0, z=0.0}
Angular-velocity vector of the wall. The value may be specified as a table of three named (x,y,z) components.

centre

Vector3, default: {x=0.0, y=0.0, z=0.0}
Axis about which the wall surface rotates. The value may be specified as a table of three named (x,y,z) components.

4.1.8. Rotating-surface, adiabatic

WallBC_RotatingSurface_Adiabatic:new{r_omega, centre, label, group}

is used where we want convective and viscous effects to impose a specified velocity condition but no heat transfer on a circular or cylindrical fixed surface. By rotating we mean that the (curved) wall is moving tangential to the block boundary such that the gas velocity is nonzero and aligned the edge of the (fixed) grid. We need to set config.viscous=true to make this boundary condition fully effective.

Twall

The fixed wall temperature, in degrees K.

r_omega

Vector3, default: {x=0.0, y=0.0, z=0.0}
Angular-velocity vector of the wall. The value may be specified as a table of three named (x,y,z) components.

centre

Vector3, default: {x=0.0, y=0.0, z=0.0}
Axis about which the wall surface rotates. The value may be specified as a table of three named (x,y,z) components.

4.2. In-flow

4.2.1. Simple supersonic

InFlowBC_Supersonic:new{flowState, x0=0.0, y0=0.0, z0=0.0, r=0.0}

is used where we want to specify a fixed (supersonic) in-flow condition that gets copied into the ghost cells each time step.

flowState

FlowState object that has been constructed earlier in your script.

Optional parameters x0, y0, z0 and r are used to set a virtual source flow. This is intended to model the flow from a conical nozzle, where the nominal flow state is at distance r from the virtual source located at x0, y0 and z0. The actual flow state on the inflow faces is computed as a perturbation to this nominal flow state.

4.2.2. Static profile

InFlowBC_StaticProfile:new{fileName, match}

is used where we want to specify an inflow condition that might vary in a complicated manner across the boundary. Data for the flow condition, on a per-cell basis, is contained in the specified file. It may be that the file is obtained from an earlier simulation, with a post-processing option like --extract-line used to write the file entries. Matching of the ghost cells to particular entries in the data file is controlled by the match string value, where the default is to match to the nearest location on all three coordinates of the ghost-cell position match="xyz-to-xyz".

Other possible values are
  • "xyA-to-xyA" For 2D or 3D simulations, don’t care about z-component of position.

  • "AyA-to-AyA" For 2D or 3D simulations, care about the y-component of position only.

  • "xy-to-xR" Starting with a profile from a 2D simulation, map it to a radial profile in a 3D simulation, considering the x-component of the position of the ghost cells.

  • "Ay-to-AR" Starting with a profile from a 2D simulation, map it to a radial profile in a 3D simulation, ignoring the x-component of the position of the ghost cells.

4.2.3. Transient

InFlowBC_Transient:new{fileName}

is used where we want to specify the time-varying inflow condition at the boundary. Data for the inflow condition, at particular time instants and assumed uniform across the full boundary, is contained in the specified file. The user needs to write this file according to the expected format encoded in the FlowHistory class, found toward the end of the flowstate.d module. Each data line will have the following space-delimited items:

time velx vely velz p T mass-fractions Tmodes (if any)

4.2.4. Constant flux

InFlowBC_ConstFlux:new{flowState, x0=0.0, y0=0.0, z0=0.0, r=0.0}

is used where we want to specify directly the fluxes of mass, momentum and energy across the boundary faces. The fluxes are computed from the supplied FlowState.

Optional parameters x0, y0, z0 and r are used to set a virtual source flow. This is intended to model the flow from a conical nozzle, where the nominal flow state is at distance r from the virtual source located at x0, y0 and z0. The actual flow state on the inflow faces is computed as a perturbation to this nominal flow state.

4.2.5. Shock-fitting

InFlowBC_ShockFitting:new{flowState, x0=0.0, y0=0.0, z0=0.0, r=0.0}

is used where we want to have the inflow boundary be the location of a bow shock. The fluxes across the boundary are computed from the supplied flow condition and the boundary velocities are set to follow the shock. Note that we need to set config.moving_grid=true and select an appropriate gas-dynamic update scheme for the moving grid.

Optional parameters x0, y0, z0 and r are used to set a virtual source flow. This is intended to model the flow from a conical nozzle, where the nominal flow state is at distance r from the virtual source located at x0, y0 and z0. The actual flow state on the inflow faces is computed as a perturbation to this nominal flow state.

4.2.6. Isentropic from stagnation

InFlowBC_FromStagnation:new{stagnationState, fileName,
  direction_type, direction_x, direction_y, direction_z,
  alpha, beta, mass_flux, relax_factor}

is used where we want a subsonic inflow with a particular stagnation pressure and temperature and a velocity direction at the boundary. Note that many of the fields are shown with their default values, so you don’t need to specify them. When applied at each time step, the average local pressure across the block boundary is used with the stagnation conditions to compute a stream-flow condition. Depending on the value for direction_type, the computed velocity’s direction can be set

  • "normal" to the local boundary,

  • "uniform" in direction and aligned with direction vector whose components are direction_x, direction_y and direction_z

  • "radial" radially-in through a cylindrical surface using flow angles alpha and beta, or

  • "axial" axially-in through a circular surface using the same flow angles.

For the case with a nonzero value specified for mass_flux, the current mass flux (per unit area) across the block face is computed and the nominal stagnation pressure is incremented such that the mass flux across the boundary relaxes toward the specified value. Note that when we select a nonzero mass flux, we no longer control the stagnation pressure. This will be adjusted to give the desired mass flux. The value for relax_factor adjusts the rate of convergence for this feedback mechanism.

Note, that for multi-temperature simulations, all of the temperatures are set to be the same as the transrotational temperature. This should usually be a reasonable physical approximation because this boundary condition is typically used to simulate inflow from a reservoir, and stagnated flow in a reservoir has ample time to equilibriate at a common temperature. The implementation of this boundary condition may not be time accurate, particularly when large waves cross the boundary, however, it tends to work well in the steady-state limit.

When mass_flux is zero and fileName is left as the default empty string, the specified FlowState is used as a constant stagnation condition. This may be modified by a user-defined function if fileName is a non-empty string that give the name of a Lua script containing a function with the name stagnationPT On every boundary condition application, this function receives a table of data (including the current simulation time) and returns values for stagnation pressure and temperature.

Here is a minimal example:

function stagnationPT(args)
   -- print("t=", args.t)
   p0 = 500.0e3 -- Pascals
   T0 = 300.0 -- Kelvin
   return p0, T0
end

The intention is that the user may program the stagnation pressure as more interesting functions of time.

4.3. Out-flow

4.3.1. Simple flux

OutFlowBC_Simple:new{}
or
OutFlowBC_SimpleFlux:new{}

is used where we want a (mostly) supersonic outflow condition. It should work with subsonic outflow as well, however, remember that you are deliberately ignoring information that may propagate into the domain from the real (physical) region that you have truncated. The outflow flux is determined from the flow state in the cell just inside the boundary. If the velocity in that cell tries to produce an influx of mass, the flux calculation switches to that of an impermeable wall.

4.3.2. Simple extrapolation

OutFlowBC_SimpleExtrapolate:new{xOrder}

is used where we want a (mostly) supersonic outflow condition. Flow data is effectively copied (xOrder=0) or linearly-extrapolated (xOrder=1) from just inside the boundary to the ghost cells just outside the boundary, every time step. In subsonic flow, this can lead to physically invalid behaviour. If you encounter strange flow behaviour that seems to start at this boundary and propagate upstream into your flow domain, try extending your simulated flow domain such that you eventually have an outflow boundary across which nothing exciting happens.

4.3.3. Fixed pressure

OutFlowBC_FixedP:new{p_outside}

is used where we want something like OutFlowBC_Simple but with a specified back pressure. This can be analogous to a vacuum pump that removes gas at the boundary to maintain a fixed pressure in the ghost cells.

4.3.4. Fixed pressure and temperature

OutFlowBC_FixedPT:new{p_outside, T_outside}

is like OutFlowBC_FixedP, above, but also sets the temperature in the ghost cells.

4.4. In-flow and Out-flow

These boundary conditions can be used on a boundary that allows inflow and outflow.

4.4.1. Ambient

InOutFlowBC_Ambient:new{flowState}

is used where we want to allow the flow to go in or out across the boundary and you are assuming an ambient flow state outside of the domain. If the velocity of the cell just inside the boundary is directed outward across the boundary face, then the interior flow properties are copied into the ghost-cells. If the cell’s flow velocity is directed inward, then the specified (ambient) flow state is copied into the ghost cells.

flowState

FlowState object that has been constructed earlier in your script.

4.4.2. DualState

InOutFlowBC_DualState:new{flowState1, flowState2, p, n}

is used where we want to allow the flow to go in or out across the boundary and you are assuming two flow states are present outside of the domain. Outside the flow domain, the flow states exist either side of a plane defined by the point p and the normal vector n. flowState1 exists on the side of the plane in the direction of n, while flowState2 is on the other side of that plane.

An example of use is in a domain that has an oblique shock processing the nominal free-stream flow. The flow before the shock would be flowState1 and after the shock flowState2, with the oblique shock being coincident with the plane.

If the velocity of the cell just inside the boundary is directed outward across the boundary face, then the interior flow properties are copied into the ghost-cells. If the cell’s flow velocity is directed inward, then relevant flow state is copied into the ghost cells. The relevant flow state is selected based on the location of a few smaple points of the boundary face. If the face sits across the plane, a blending of the flow states is used.

flowState1

FlowState object that has been constructed earlier in your script. On side of plane in the direction of n.

flowState2

FlowState object that has been constructed earlier in your script. On side of plane in opposite to the direction of n.

p

Vector3, default: {x=0.0, y=0.0, z=0.0}
A point on the plane that separates the flow states. The value may be specified as a table of three named (x,y,z) components.

n

Vector3, default: {x=0.0, y=1.0, z=0.0}
Normal vector to the plane. The value may be specified as a table of three named (x,y,z) components.

4.5. Inter-block Exchange

4.5.1. Full block-face

ExchangeBC_FullFace:new{otherBlock, otherFace, orientation,
  reorient_vector_quantities, Rmatrix}

Usually, this boundary condition is applied implicitly, by calling the function identifyBlockConnections, for cases where one structured-grid block interfaces with another and the block boundaries are cleanly aligned, however, it can be applied manually for cases where you want the flow to be plumbed from one block face into another and the blocks are not geometrically aligned. A non-unity transformation matrix, Rmatrix, can be provided for cases where the flow vector quantities need to be reoriented when they are copied from the other boundary to this one.

Note that this boundary condition is only for structured-grid blocks. If one or both of the blocks to be joined is based on an unstructured-grid, you will need to use the following MappedCell flavour of the exchange boundary condition. Also, if you have two structured-grid faces but the cells along each of the faces do not align, you will also need to use the MappedCell flavour.

otherBlock

FluidBlock object, default: nil
A reference to the other block that is joined to this boundary.

otherFace

FaceLabel, default: nil
An enum value specifying the face on the other block, from which cell data will be copied. Possible values are north, east, south, west, bottom, top.

orientation

int: default: -1
Although there is only one possible orientation for face-to-face connections in 2D, there are 4 possible ways to (rotationally) orient a face-to-face connection in 3D.

reorient_vector_quantities

boolean, default: false
If true, vector quantities are multiplied by the Rmatrix as they are copied from the source cell into the destination ghost cell.

Rmatrix

array of float, default: {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}
This is a 3-by-3 matrix written in row-major format. It is applied, in the usual matrix-vector multiply manner, to vector quantities such as velocity and magnetic field. The user is responsible for computing appropriate coefficients.

4.5.2. Mapped cell

ExchangeBC_MappedCell:new{transform_position, c0, n, alpha, delta,
  list_mapped_cells, reorient_vector_quantities, Rmatrix}

is something like the ExchangeBC_FullFace boundary condition but with a mapping of destination(ghost)-cell location to source-cell location. It allows us to stitch boundaries together, even if the cells do not align, one-for-one. The position of the source cell is computed by taking the position of the ghost cell, computing the solid-body rotation of alpha radians about the axis n through the point c0, then adding a displacement delta. This will accommodate general rigid-body transformations.

c0

Vector3, default: {x=0.0, y=0.0, z=0.0}
Centre of rotation for position transformation.

n

Vector3 default: {x=0.0, y=0.0, z=1.0}
Axis of rotation for the position transformation.

alpha

float, default: 0.0
Angle of rotation for the position transformation. Right-hand rule gives the sense of direction.

delta

Vector3, default: {x=0.0, y=0.0, z=0.0}
Translational displacement for the position transformation.

list_mapped_cells

boolean, default: false
Flag to indicate whether we want the program to list the indices of the mapped cells as they are located at run time.

reorient_vector_quantities

boolean, default: false
If true, vector quantities are multiplied by the Rmatrix as they are copied from the source cell into the destination ghost cell.

Rmatrix

array of float, default: {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}
This is a 3-by-3 matrix written in row-major format. It is applied, in the usual matrix-vector multiply manner, to vector quantities such as velocity and magnetic field. The user is responsible for computing appropriate coefficients.

4.6. User Defined

This is a get-out-of-jail boundary condition that allows you to do anything you wish to code (in Lua). Please read the Appendix on these boundary conditions in the PDF User Guide before your first attempt to use them.

4.6.1. Setting ghost cells

UserDefinedGhostCellBC:new{fileName, label, group}

is used to allow the user to define the ghost-cell flow properties and interface properties at run time. This is done via a set of functions defined by the user, written in the Lua programming language, and provided in the specified file.

4.6.2. Setting fluxes

UserDefinedFluxBC:new{fileName, funcName, label, group}

is used to allow the user to define the interface convective-fluxes at run time. This is done via a function defined by the user, written in the Lua programming language, and provided in the specified file. If the user does not specify the function name, convectiveFlux is used as the default name.

4.6.3. Full block-face followed by user-defined

ExchangeBC_FullFacePlusUDF:new{otherBlock, otherFace, orientation,
  reorient_vector_quantities, Rmatrix, fileName, label, group}

is used when you might conditionally want to exchange block-boundary data or do something else altogether. This boundary condition allows that by first doing a FullFace exchange of data and then calling upon your user-defined functions (as for UserDefinedGhostCellBC) to conditionally overwrite the data. This turns out to be a convenient way to implement diaphragm models for shock-tunnel simulations. Note that this boundary condition can work across MPI tasks but is only implemented for structured-grid blocks.

5. FlowSolution

In a preparation script, the FlowSolution object is used to initialise the domain with a flow field from another (completed) job.

FlowSolution:new{jobName, dir, tindx, nBlocks, gindx, make_kdtree}
jobName

string, no default
base file name for completed simulation

dir

string, no default
directory containing completed simulation
It is often easy to use relative paths if the simulation is nearby.

tindx

int or string, no default
time index at which to read flow field from completed simulation
provide an integer for the time index, or to get the last flow field, one can supply the special string "last".

nBlocks

int, no default
number of blocks in completed simulation.
Note that if you provide a value smaller than the total number of blocks in the completed simulation, then only blocks up to nBlocks-1 will be used to initialise the flow field. A use case for this would if one only want to use some inflow portion of the domain when intialising a new simulation.

gindx

int or string, default = -1
grid time index at which to read the grid associated with the flow field
For the default value of -1, the grid time index is synced with the tindx choice. This is likely what you want for moving grids, and for non-moving grids this will just pick up the grid at grid time index 0.
A special use for this gindx parameter is when attempting to initialise a moving grid simulation from a non-moving grid simulation. You can set the grid time index explicitly to 0 so that the non-moving flow solution is read in with its initial grid at grid time index 0.

make_kdtree

int, default = 0
when make_kdtree is selected (=1) the search speed for locating cells is considerably improved. It does this by using a kd tree filtering and sorting approach. There is, however, a startup cost for building the kd tree.

Example of use:

initial = FlowSolution:new{jobName='sphere', dir='../inviscid',
                           tindx=10, nBlocks=4, make_kdtree=1}

6. Geometric elements

The basic geometric elements available are points (Vector3), paths, surfaces and volumes.

6.1. Points

In the Eilmer input, any point can be specified as a Lua table with x, y and z components or as a Vector3 object. If you are declaring points with no further need to perform vector operations, then the table format is preferred (to reduce the syntactic noise in your input scripts). To demonstrate these alternatives, here are examples of two valid inputs for constructing a line from (1,2) to (4,5):

ab = Line:new{p0={x=1.0, y=2.0}, p1={x=4.0, y=5.0}}
ab2 = Line:new{p0=Vector3:new{x=1.0, y=2.0}, p1=Vector3:new{x=4.0, y=5.0}}

When using a Lua table, any missing x, y or z components will default to a value of 0.0.

Obviously, the form with points as Lua tables is more compact. Why then would one want to use the long-form Vector3 objects? Well, these are useful when you need to do calculations involving vector arithmetic. This vector arithmetic (as vector operations) has been built for you in the Eilmer ecosystem. Read on in the next section.

6.2. Vector3

This class defines geometric vector objects with three Cartesian components: x, y and z. The Vector3 constructor accepts these components in a table, or another Vector3 object.

p0 = Vector3:new{x=1.0, y=2.0, z=3.0}
p1 = Vector3:new{p0}

6.2.1. Vector3 expressions

A number of methods have been defined so that you can write arithmetic expressions that involve Vector3 objects.

Sample expressions available for Vector3 objects include:

p2 = -p0           -- negative copy --> Vector3:new{x=-1.0, y=-2.0, z=-3.0}
p2 = p0 + p1       -- addition --> Vector3:new{x=2.0, y=4.0, z=6.0}
p2 = p0 - p1       -- subtraction --> Vector3:new{x=0.0, y=0.0, z=0.0}
p2 = p0 * 3.0      -- scaling --> Vector3:new{x=3.0, y=6.0, z=9.0}
p2 = 3.0 * p0      -- scaling --> Vector3:new{x=3.0, y=6.0, z=9.0}
p2 = p0 / 3.0      -- scaling --> Vector3:new{x=1./3, y=2./3, z=3./3}
p2:normalize()     -- scales to unit magnitude -->
                   --     Vector3:new{x=0.267261, y=0.534522, z=0.801783}
p3 = unit(p0)      -- new unit vector of p0 -->
                   --     Vector3:new{x=0.267261, y=0.534522, z=0.801783}
a = vabs(p3)       -- magnitude --> 1.0
p1 = Vector3:new{p0}
b = dot(p0, p1)    -- dot product --> 14.0
p3 = cross(p0, p1) -- cross product --> Vector3:new{x=0.0, y=0.0, z=0.0}

6.3. Path elements

Path objects are typically used to define boundaries for surfaces. A Path object can be called to evaluate a point at parameter t, where the parametric range is 0.0 to 1.0.

The constructors for Path objects include:

Line:new{p0, p1}

Defines a straight line from point p0 (t=0) to point p1 (t=1).

Bezier:new{points}

Defines a Bézier curve from the sequence of points.

NURBS:new{points, weights, knots, degree}

Defines a NURBS from control points with weights, knot vector and degree.

Arc:new{p0, p1, centre}

Defines a circular arc from point p0 to point p0 about centre.

Arc3:new{p0, pmid, p1}

Defines a circular through three points starting at p0 through pmid ending at p1.

Polyline:new{segments}

Builds a single path from a sequence of Path objects defined in segments.

Spline:new{points}

Builds a spline of Bezier segments through the sequence of points.

Spline2:new{filename}

Builds a spline from data file specifed in filename key with points listed line by line.

LuaFnPath:new{luaFnFname}

Builds a path based on a user-supplied function (as a function of t) specified in luaFnName key.

ArcLengthParameterizedPath:new{underlying_path}

Derives path from underlying_path that has a uniformly-distributed set of points with parameter t.

6.4. Surfaces (Patches)

Surfaces or patches are elements that are parameterised in with two coordinates. In our definition, those coordinates are r and s, and the range for each is from 0.0 to 1.0. The surface itself can be two-dimensional or three-dimensional. Two-dimensional surfaces are often used in structured grid generation for 2D simulations. For this use case, we often refer to these entities are patches. Three-dimensional surfaces can be used to represent faces on 3D volumes, so they often turn up in 3D grid generation.

In the following, the syntax for various surface types is given and described. We have also produced some grids from the various surface types to give you some idea of the effect of surface type choice. Three example geometries, covering a modest range of use cases, have been chosen for the demonstration:

  1. a blunt-body grid

  2. a converging-diverging nozzle grid

  3. a duct-like grid

The examples appear at the end of this section.

6.4.1. Coons patch

There are two constructors available for a Coons patch. One accepts four Path objects, the other accepts four corner points ( Vector3 s).

CoonsPatch:new{north, south, east, west}
CoonsPatch:new{p00, p10, p11, p01}

6.4.2. Area-orthogonal patch

The area-orthogonal patch attempts to equalise the areas of cells and maintain orthogonality at the edges. The constructor is:

AOPatch:new{north, south, east, west, nx=10, ny=10}

The nx and ny parameters control the background grid used to build an area-orthogonal patch. The default is 10 background-cells in each direction. If you have boundaries that have sections high curvature you may need to increase the number of background cells.

Note that with default settings the blunt-body grid results in negative cells near the shoulder and the nozzle grid has cut the corner of the boundary. This behaviour can be improved if some care is take in selecting the background grid. The next example, for the nozzle, increased the background grid in the x-direction to nx=40 to give a better representation of the change in geometry in the streamwise direction. This is show here.

6.4.3. Channel patch

The channel patch constructor is:

ChannelPatch:new{north, south, ruled=false, pure2D=false}

The channel patch builds an interpolated surface between two bounding paths. Those bounding paths are specified as north and south. Bridging paths are built between the bounding paths. By default, the cubic Bézier curves are used as the bridging paths. This choice provides good orthoganality of grid lines departing the bounding paths.

If ruled=true is set, then the bridging curves are straight-line segments.

If pure2D=true is set, then all z-ordinates are set to 0.0. This might be useful for purely two-dimensional constructions.

6.4.4. Ruled surface

The constructor for a ruled surface is:

RuledSurface:new{edge0, edge1, ruled_direction='s', pured2D=false}

The ruled surface is used to build a surface connecting to opposing paths with straight-line segments. For example, a ruled surface based on south path (sPath) and a north path (nPath) could be constructed as:

mySurf1 = RuledSurface:new{edge0=sPath, edge1=nPath}

In this case, the default ruled_direction='s' is appropriate.

For using west and east paths as the edge paths, the constructor is:

mySurf2 = RuledSurface:new{edge0=wPath, edge1=ePath, ruled_direction='r'}

The constructor also accepts i and j as direction designators in place of r and s.

If pure2D=true is set, then all z-ordinates are set to 0.0. This might be useful for purely two-dimensional constructions.

6.4.5. Control point patch

The control point patch is an algebraic surface type with interior control of grid lines by means of control points. There are two constructors available. One accepts the control points directly, as created independently by the caller. The other constructor acceps the number of control points in each logical dimension and a guide patch. The guide patch is used to create and distribute the control points. The constructors are:

ControlPointPatch:new{north, south, east, west, control_points}
ControlPointPatch:new{north, south, east, west,
                      ncpi, ncpi, guide_patch}

For the first constructor, the control_points are listed as a table of tables such that the points are given as a 2D array.

For the second constructor, ncpi and ncpj are integers specifying the numbers of control points in the i and j dimensions, respectively. The guide_patch is given as a string; supported options are:

  • "coons"

  • "aopatch"

  • "channel"

  • "channel-e2w"

The channel patch is defined assuming that the north and south boundaries are the bounding walls of the channel. When using a channel patch as a guide patch for the ControlPointPatch, it is sometimes convenient to have the channel oriented with east and west boundaries as bounding walls. The guide_patch set to "channel-e2w" gives that option.

More information about the control point patch is available in the technical note Control Point Form Grid Generation.

6.4.6. Patch examples for 2D grid generation

These examples show how the various patch types work on a common set of bounding geometries.

A blunt body grid (FireII)
fire2 coons patch grid
fire2 AO patch grid
fire2 channel patch grid

(a) Coons patch

(b) Area-orthogonal patch

(c) Channel patch

fire2 ruled surface grid
fire2 ctrl pt patch grid
fire2 ctrl pt patch grid p pts

(d) Ruled surface

(e) Control point patch

(f) Control point patch with control points

Note that with default settings the blunt-body grid generated with the Area-orthogonal patch results in negative cells near the shoulder. To remedy this, increase the ny value for number of points in the background grid to get better resolution of the geometry change around the shoulder.

A nozzle grid
nozzle coons patch grid
nozzle AO patch grid

(a) Coons patch

(b) Area-orthogonal patch

nozzle channel patch grid
nozzle ruled surface grid

(c) Channel patch

(d) Ruled surface

nozzle ctrl pt patch grid
nozzle ctrl pt patch grid p pts

(e) Control point patch

(f) Control point patch with control points

Here the Area-orthogonal patch with default settings as cut the corner of the boundary on the nozzle domain. This behaviour can be improved if some care is take in selecting the background grid. We show this here by increasing the background grid in the x-direction to nx=40 to give a better representation of the change in geometry in the streamwise direction.

nozzle AO patch grid 2
A duct grid

This duct grid mimics the example in Eiseman (1988) A control point form of algebraic grid generation. The channel patch and ruled surface are not appropriate for this domain because the curved boundaries on both the north and east edges. As such, they are not generated for this example.

duct coons patch grid
duct AO patch grid

(a) Coons patch

(b) Area-orthogonal patch

duct ctrl pt patch grid
duct ctrl pt patch grid p pts

(c) Control point patch

(d) Control point patch with control points

7. Grid objects, methods and functions

7.1. StructuredGrid objects, methods and functions

Here is the list of methods and functions that are specific to structured grids. The more general set of methods and functions that apply to all grids (structured and unstructured) are documented in Grid objects, methods and functions.

StructuredGrid:joinGrid(otherGrid, joinLocation)

Use this method to join another grid onto the current object. If successful, the current object is modified at the end of the operation. It will now contain both grids.

The join grid operation assumes that the edge/face for joining is full-face matching and that the grid points are C0 continuous. The operation does not do extensive checks for badly formed joins. It relies on the caller to use this in a sensible manner.

otherGrid

second StructuredGrid to be joined

joinLocation

a string specifying at which end of the current grid the other grid should be joined
allowable options are "east", "imax", "north" and "jmax". Note that the first two are equivalent, and the last two are equivalent.

8. Running a simulation

The overall process of simulation of a gas flow is divided into three stages:

(1) preparation of the defining data,
(2) running the main simulation program, and
(3) postprocessing of the output data from the simulation.

Information is passed from one stage to the next in data files.

8.1. Pre-processing

To start a simulation, we need a gas model, maybe a chemical-kinetics model for reacting flow, simulation configuration and control files, and the initial grid and flow state. Generation of all of these data is the activity of the pre-processing (or preparation) stage.

8.1.1. prep-gas : gas model preparation

The prep-gas program is used to take a brief description of the desired gas model used in the flow simulation and produce a detailed configuration gas model configuration file for use by Eilmer at pre-processing and simulation stages. Its usage is shown here. Generally, one uses prep-gas in the first mode shown: with two arguments. The second mode simply lists available species in the Eilmer database and exits.

Usage:
 > prep-gas input output

   input    : a gas input file with user selections
   output   : detailed gas file in format ready for Eilmer4.

 > prep-gas --list-available-species

8.1.2. prep-chem : chemistry scheme preparation

prep-chem is used to take user created description of a chemistry scheme written in Lua and generate a detailed configuration file for eilmer to use at run-time. The use of prep-chem is shown here.

Usage:
 > prep-chem [--compact] gmodelfile cheminput output

   gmodelfile  : a gas model file is required as input for context
   cheminput   : input chemistry file in Lua format.
   output      : output file in format ready for Eilmer4.

Options:
   --compact   : produce a text file called 'chem-compact-notation.inp'
                 which is used to configure a GPU chemistry kernel.

8.1.3. e4shared --prep : pre-processing for flow simulation

Once the gas and chemistry files are prepared, the flow simulation process in divided into three stages, with preparation of the configuration, grid and initial flow-state files being the first.

 > e4shared --prep --job=name [--verbosity=<int>]

will start with the user-supplied description of the simulation in the file name.lua and produce the set of data files that are needed to run a simulation. The integer value for verbosity defaults to 0, for minimum commentary as the preparation calculations are done. If things are going badly, setting a higher value may give more information during the calculations, however, it is more likely that increased verbosity is more useful during the simulation (run-time) stage.

8.1.4. e4loadbalance : load-balancing when n_blocks exceeds n_processors

For parallel simulations with MPI, it might be useful to load balance the work by stacking several computational blocks onto a single MPI process. This load balance tool will estimate load based on number of cells in a block and try to distribute the work in an even manner. The user should decide how many processes are available. A typical use is when the number of blocks significantly exceeds the number of processes available. As an example, a simulation of 100 computational blocks could be distributed onto a 32-core workstation using this load balance tool. The command would be:

 > e4loadbalance --job=name --ntasks=32

This works for multi-node simulations also.

To inspect the result of the load balance, you can look at the mpimap file that is written to config directory: config/<name>.mpimap. This file will show the mapping of blocks on to process IDs. (It could be hand-edited if you know what you’re doing.)

8.2. Run-time

Now that the initial flow and grid files exist, it is time to run the main transient flow simulation process. The transient-flow simulation program comes in two main flavours. e4shared is for using one or more cores on your workstation for calculation using shared-memory parallelism. e4mpi is for running a distributed-memory calculation using the Message-Passing Interface (MPI). When getting started with small-scale simulations, e4shared is the simpler program to use. Once you have graduated to running large, many-block simulations, it may be good to start using e4mpi across the nodes of a cluster computer.

8.2.1. Run a simulation, shared-memory

 > e4shared --run --job=name [OPTIONS]

with the relevant options being

  --tindx-start=<int>|last|9999      defaults to 0
  --next-loads-indx=<int>            defaults to (final index + 1) of lines
                                     found in the loads.times file
  --max-cpus=<int>                   (e4shared) defaults to 8 on this machine
  --max-wall-clock=<int>             in seconds
  --report-residuals                 include residuals in console output

8.2.2. Run a simulation, MPI

 > mpirun -np <ntask> e4mpi --run --job=name [OPTIONS]

with the relevant options being

  --tindx-start=<int>|last|9999      defaults to 0
  --next-loads-indx=<int>            defaults to (final index + 1) of lines
                                     found in the loads.times file
  --threads-per-mpi-task=<int>       (e4mpi) defaults to 1
  --max-wall-clock=<int>             in seconds
  --report-residuals                 include residuals in console output

8.3. Post-processing

Once you have a large collection of numbers defining the flow field at various instances in its history, it will be stored in files with a format specific to Eilmer and mostly unknown to many data visualization tools. To make your data more accessible, you may rewrite it into formats known by the visualization tools or you may slice it into various subsets.

8.3.1. Post-processing using predefined actions

 > e4shared --post --job=name [OPTIONS]

with the relevant options being

  --list-info                        report some details of this simulation
  --tindx-plot=<int>|all|last|9999   defaults to last
  --add-vars="mach,pitot"            add variables to the flow solution data
                                     (just for postprocessing)
                                     Other variables include:
                                       total-h, total-p, total-T,
                                       enthalpy, entropy, molef, conc,
                                       Tvib (for some gas models)
                                       nrf (non-rotating-frame velocities)
                                       cyl (cylindrical coordinates: r, theta)
  --ref-soln=<filename>              Lua file for reference solution
  --vtk-xml                          produce XML VTK-format plot files
  --binary-format                    use binary within the VTK-XML
  --tecplot                          write a binary szplt file for Tecplot
  --tecplot-ascii                    write an ASCII (text) file for Tecplot
  --plot-dir=<string>                defaults to plot
  --output-file=<string>             defaults to stdout
  --slice-list="blk-range,i-range,j-range,k-range;..."
                                     output one or more slices across
                                     a structured-grid solution
  --surface-list="blk,surface-id;..."
                                     output one or more surfaces as subgrids
  --extract-streamline="x,y,z;..."   streamline locus points
  --track-wave="x,y,z(,nx,ny,nz);..."
                                     track wave from given point
                                     in given plane, default is n=(0,0,1)
  --extract-line="x0,y0,z0,x1,y1,z1,n;..."
                                     sample along a line in fluid domain
  --extract-solid-line="x0,y0,z0,x1,y1,z1,n;..."
                                     sample along a line in solid domain
  --compute-loads-on-group=""        group tag
  --probe="x,y,z;..."                locations to sample flow data
  --output-format=<string>           gnuplot|pretty
  --norms="varName,varName,..."      report L1,L2,Linf norms
  --region="x0,y0,z0,x1,y1,z1"       limit norms calculation to a box

8.3.2. Post-processing using a user-supplied script

When none of the predefined post-processing operations are suitable, you may define your own, in Lua. The Eilmer4 program provides a number of service functions to the Lua interpreter for loading grid and flow files and accessing the data within the loaded grids and flow blocks. This is probably the least-well-defined activity associated with a simulation, so an interest in experimentation could be rewarding.

 > e4shared --custom-post --script-file=name.lua

9. Restarting from a snapshot

You may instruct Eilmer to make snapshots of the solution at every N steps. You can also choose how many snapshots to keep as the simulation progresses. This facility is useful for large simulations that are running on cluster computers. If something happens during the job run, having a snapshot might give you a point to return to such that you haven’t lost too much compute time. To configure snapshots in your input script, add the following lines:

config.snapshot_count = 1000
config.number_total_snapshots = 2

In this example, snapshots are written every 1000 steps and a total of two snapshots are kept. Say you simulation has run for 5692 steps. We would expect two snapshots in the flow/ directory:

snapshot-0000, recorded at step 4000
snapshot-0001, recorded at step 5000

When the simulation gets to 6013 steps, the snapshots in the flow directory are:

snapshot-0000, recorded at step 5000
snapshot-0001, recorded at step 6000

You can see the pattern here. The largest snapshot index corresponds to the most recent snapshot.

If your calculation is interrupted, you can restart from one of the snapshots. There is some file house-keeping required to insert the snapshot into the regular time series. The e4-prep-restart program takes care of that house-keeping effort. There are two modes for the e4-prep-restart program:

  1. to append a snapshot to the time series; or

  2. to replace a solution in the time series with a snapshot.

The first mode is more common and so is the default action. Let’s look at an example of appending one of our snapshots to the time series.

> e4-prep-restart --job=my_job --snapshot=1

This will take the snapshot at index 1 and append it to the time series. You will get a message to the screen informing you of the new time index in the series. Take note of this new index as you will likely want to use that to instruct an Eilmer restart. For example, if the new time index is 9, we could restart an Eilmer job using:

> e4shared --job=my_job --run --tindx-start=9

The other mode of operation for e4-prep-restart is to replace a solution in the time series with a snapshot. In this example, we replace the solution at index 5 with the snapshot 0:

> e4-prep-restart --job=my_job --snapshot=0 --replace=5

Since we’ve replaced index 5 in the time series, it’s likely we’d want to restart Eilmer from that point. Here is an example with MPI (on a cluster where number of processes is taken care of):

> mpirun e4mpi --job=my_job --run --tindx-start=5

The e4-prep-restart program accepts a --verbose flag that will show you the house-keeping steps being performed. Some command help is available with the --help flag.

10. Run-time utilities

10.1. Monitoring a simulation

While a calculation runs, its progress can be monitored in a couple of ways. By default, the simulation program occasionally writes to the console the time-step count, the current simulation time, the current size of time step, the current wall-clock time (WC), and estimates of the wall-clock time until the simulation reaches either the final simulation time (WCtFT) or the maximum allowed number of steps (WCtMS). For a simple simulation of a transient flow, this is often sufficent.

10.1.1. Viewing residuals

For a simulation that is intended to progress toward a steady solution, you can monitor the residuals in the flow-field by specifying the --report-residuals command-line option when starting the simulation. With this option, the residuals for mass, momentum in each of the coordinate directions, energy, the L2 norm of the collection of residuals and the mass-balance, are written to the file config/<job>-residuals.txt. This file can be monitored by writing it to another console with tail -f or its content can be rendered graphically with GNUPlot. See dgd/examples/eilmer/2D/lat-plate-hakkinen-swbli/sg-lts/ for a sample GNUPlot script for rendering of this data. The essential lines from that script are shown here.

job = "swbli"
period_seconds = 2.0
residuals_file = "./config/" . job ."-residuals.txt"
window_title = "Eilmer4 residuals for ".job
set term qt title window_title
set logscale y
bind "x" "end=1"
bind "q" "end=1"
end = 0
while( end==0 ) {
  set multiplot layout 2,2
  plot residuals_file using "step":"x-mom" with lines, \
    residuals_file using "step":"y-mom" with lines
  plot residuals_file using "step":"energy" with lines
  plot residuals_file using "step":"mass" with lines
  plot residuals_file using "step":"L2" with lines
  pause period_seconds
}

10.1.2. (Lack of) Simulation progress

For large MPI jobs, there is a risk that the part of the calculation will fail and leave other MPI tasks consuming compute resources while waiting for communication that will never come. We provide the e4monitor program that can periodically check the progress of a simulation and, if progress is not being made, terminate the calculation process forcefully.

Progress is indicated by the simulation program writing the current time step number to the file config/<job>-progress.txt. This file is overwritten each time step. On completion of a simulation, the word done is written to that file.

To use the e4monitor program in your batch scripts, start it just before using mpirun to start your simulation.

$ e4monitor --job=<jobName> ?--startup=<int>? ?--period=<int>? &

The trailing ampersand indicates that the e4monitor program has been started in the background. The default values for startup and period are 30 and 60 seconds respectively. With these values, the monitor program does not look at the progress file until 30 seconds has elapsed. It subsequently checks every 60 seconds that the time-step count has increased. Set these parameters to suit your particular simulation. If the monitor program does not see progress for a few cycles, it issues a kill command for the mpirun process. If the simulation finishes successfully, the monitor program will see the word done and will also terminate. This should give a clean exit for your batch job.

10.2. Compacting a solution

When running a transient-flow simulations with many blocks and saving the flow data at many time instants, you can end up with thousands of files. Although not usually a problem on your personal workstation, storing many small files on a shared supercomputer file system is likely to lead to trouble with your storage allocation. To mitigate this problem, there is a utility program to compact and restore job directories.

$ e4compact --job=<jobName> ?--verbose? ?--force? ?--restore?

This program delegates most of the work to the tar command and will collect the many solution and grid files into a small set of archive files. Files that are successfully archived are removed from the directories and, if the directories are then empty, they are also removed.

Notes
  • Existing archives will only be updated if the --force|-f option is specified.

  • Option --restore|-r will extract files from archives and restore the directories.

  • Option --verbose|-v will echo the tar commands as they are issued.

11. Extract Boundary Loads

Sometimes the primary objective for a fluid simulation is to calculate the loads on a particular boundary of the domain, for example, calculating the drag or heat load along the surface of a wing or reentry capsule. Eilmer allows users to easily access the boundary data on specified boundary conditions for post-processing. It is important to note that the configuration settings presented below must be set before a simulation, extracting the boundary loads data is a simulation run-time procedure.

config.write_loads

boolean, default: false
To extract boundary loads, set config.write_loads=true.

config.boundary_groups_for_loads

string, default: 'loads'
The config.boundary_groups_for_loads setting is used to select which boundaries will have their loads data written to disk. This string should match the groups string set in the corresponding boundary condition. This can be a list of strings separated by a comma, e.g. config.boundary_groups_for_loads='group1,group2'.

config.dt_loads

float, default: 1.0e-03
The boundary loads data will be written to disk repeatedly, each time this increment of simulation time has elapsed.