case.md

@page case Case Files

Case Files

Example Python case files, also referred to as input files, can be found in the examples/ directory. They print a Python dictionary containing input parameters for MFC. Their contents, and a guide to filling them out, are documented in the user manual. A commented, tutorial script can also be found in examples/3d_sphbubcollapse/.

Basic Skeleton

The skeleton for an input file may look like the following:

#!/usr/bin/env python3

import json

# Configuring case dictionary
print(json.dumps({
  # Insert case parameters here
  ...
}))

Thus, you can run your case file with Python to view the computed case dictionary that will be processed by MFC when you run:

python3 my_case_file.py

This is particularly useful when computations are done in Python to generate the case.

(Optional) Accepting command line arguments

Input files can accept command line arguments, forwarded by mfc.sh run. Consider this example from the scaling case:

import json, argparse

parser = argparse.ArgumentParser(
    prog="scaling",
    description="Weak- and strong-scaling benchmark case.",
    formatter_class=argparse.ArgumentDefaultsHelpFormatter)

parser.add_argument("--mfc", type=json.loads, default='{}', metavar="DICT",
                    help="MFC's toolchain's internal state.")
parser.add_argument("-s", "--scaling", type=str, metavar="SCALING", choices=["weak", "strong"],
                    help="Whether weak- or strong-scaling is being exercised.")

# Your parsed arguments are here
args = parser.parse_args()

The --mfc argument is a JSON string representing mfc.sh run's internal state, passed in when MFC runs your input file. It contains all the runtime information you might want from the build/run system. You can add as many additional arguments and options as you may need.

To run such a case, use the following format:

./mfc.sh run <path/to/case.py> <mfc.sh run arguments> -- <case arguments>

For example, to run the scaling case in "weak-scaling" mode:

./mfc.sh run examples/scaling/benchmark.py -t pre_process -j 8 -- --scaling weak

Parameters

Feature Compatibility

Before diving into parameter details, check the @ref case_constraints "Feature Compatibility Guide" to understand:

• Which features work together (MHD, bubbles, phase change, etc.)
• Common configuration patterns with copy-paste examples
• Requirements for each model equation and Riemann solver

Parameter Lookup:

• CLI search: ./mfc.sh params <query> - Search ~3,300 parameters from the command line
• Full reference: @ref parameters "Case Parameters" - Complete parameter documentation

There are multiple sets of parameters that must be specified in the python input file:

1. Runtime Parameters
2. Computational Domain Parameters
3. Patch Parameters
4. Immersed Boundary Patches
5. Fluid Material's Parameters
6. Simulation Algorithm Parameters
7. Formatted Database and Structure Parameters
8. (Optional) Acoustic Source Parameters
9. (Optional) Ensemble-Averaged Bubble Model Parameters
10. (Optional) Velocity Field Setup Parameters
11. (Optional) Phase Change Parameters
12. (Optional) Artificial Mach Number Parameters

Items 8, 9, 10, 11 and 12 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, initial velocity field setup, phase change, artificial Mach number respectively. Definition of the parameters is described in the following subsections.

Enumerated parameters accept named values as well as integer codes: "riemann_solver": "hllc" is equivalent to "riemann_solver": 2. Defined names appear in each parameter's table entry (e.g. 1 (hll), 2 (hllc)). Existing case files can be rewritten to named syntax with ./mfc.sh validate case.py --migrate.

1. Runtime {#sec-runtime}

Parameter	Type	Description
run_time_info	Logical	Output run-time information
rdma_mpi	Logical	(GPUs) Enable RDMA for MPI communication.
active_box	Logical	Enable causal-envelope active-box restriction of the RHS compute window.
case_dir	String	Case directory path
old_grid	Logical	Use grid from previous simulation
old_ic	Logical	Use initial conditions from previous simulation
t_step_old	Integer	Time step to restart from
n_start_old	Integer	Starting index from previous simulation

• run_time_info generates a text file that includes run-time information including the CFL number(s) at each time-step.
• active_box enables the causal-envelope active-box optimization, restricting the RHS compute window to the region where the solution deviates from a uniform ambient state. Requires WENO reconstruction (recon_type = 1) and SSP-RK3 time stepping (time_stepper = 3). Incompatible with immersed boundaries, acoustic sources, body forces, Lagrangian bubbles, phase change, and the IGR solver.
• rdma_mpi optimizes data transfers between GPUs using Remote Direct Memory Access (RDMA). The underlying MPI implementation and communication infrastructure must support this feature, detecting GPU pointers and performing RDMA accordingly.

2. Computational Domain {#sec-computational-domain}

Parameter	Type	Description
x[y,z]_domain%%beg[end]	Real	Beginning [ending] of the $x$[y,z]-direction domain
stretch_x[y,z]	Logical	Stretching of the mesh in the $x$[y,z]-direction
a_x[y,z]	Real	Rate at which the grid is stretched in the $x$[y,z]-direction
x[y,z]_a	Real	Beginning of the stretching in the negative $x$[y,z]-direction
x[y,z]_b	Real	Beginning of the stretching in the positive $x$[y,z]-direction
loops_x[y,z]	Integer	Number of times to recursively apply grid stretching
cyl_coord	Logical	Cylindrical coordinates (2D: Axisymmetric, 3D: Cylindrical)
m	Integer	Number of grid cells in the $x$-coordinate direction
n	Integer	Number of grid cells in the $y$-coordinate direction
p	Integer	Number of grid cells in the $z$-coordinate direction
pref	Real	Reference pressure
rhoref	Real	Reference density

The parameters define the boundaries of the spatial and temporal domains, and their discretization that are used in simulation.

• [x,y,z]_domain%[beg,end] define the spatial domain in $x$, $y$, and $z$ Cartesian coordinates:

\f[ x \in \left[ x_domain%beg, ; x_domain%end \right], \quad y \in \left[ y_domain%beg, ; y_domain%end \right], \quad z \in \left[ z_domain%beg, ; z_domain%end \right] \f]

• $m$, $n$, and $p$ define the number of finite volume cells that uniformly discretize the domain along the $x$, $y$, and $z$ axes, respectively. Note that the actual number of cells in each coordinate axis is given as $[m,n,p]+1$. For example, $m=n=p=499$ discretizes the domain into $500^3$ cells. When the simulation is 2D/axi-symmetric or 1D, it requires that $p=0$ or $p=n=0$, respectively.

• stretch_[x,y,z] activates grid stretching in the $[x,y,z]$ directions. The grid is gradually stretched such that the domain boundaries are pushed away from the origin along a specified axis.

• a_[x,y,z], [x,y,z]_a, and [x,y,z]_b are parameters that define the grid stretching function. When grid stretching along the $x$ axis is considered, the stretching function is given as:

\f[ x_{cb,stretch} = x_{cb} + \frac{x_{cb}}{a_x} \Bigg[ \mathrm{log}\left[\mathrm{cosh} \left( \frac{a_x(x_{cb}-x_a)}{L} \right) \right] + \mathrm{log}\left[\mathrm{cosh} \left( \frac{a_x(x_{cb}-x_b)}{L} \right) \right] -2 \mathrm{log}\left[\mathrm{cosh} \left( \frac{a_x(x_b-x_a)}{2L} \right) \right] \Bigg] \f]

where x_cb and x_[cb,stretch] are the coordinates of a cell boundary at the original and stretched domains, respectively. L is the domain length along the x axis: L=x_domain%%end-x_domain%%beg. Crudely speaking, x_a and x_b define the coordinates at which the grid begins to get stretched in the negative and positive directions along the $x$ axis, respectively. $a_x$ defines the smoothness of the stretching. Stretching along the $y$ and $z$ axes follows the same logistics. Optimal choice of the parameters for grid stretching is case-dependent and left to the user. loops_x[y,z] defines the number of times the grid stretching function is applied and has a default value of one.

• cyl_coord activates cylindrical coordinates. The domain is defined in $x$-$y$-$z$ cylindrical coordinates, instead of Cartesian coordinates. Domain discretization is accordingly conducted along the axes of cylindrical coordinates. When $p=0$, the domain is defined on $x$-$y$ axi-symmetric coordinates. In both Coordinates, mesh stretching can be defined along the $x$- and $y$-axes. MPI topology is automatically optimized to maximize the parallel efficiency for given choice of coordinate systems.

3. Patches {#sec-patches}

Parameter	Type	Analytical Definition	Description
num_patches	Integer	Not Supported	Number of initial condition geometric patches.
num_fluids	Integer	Not Supported	Number of fluids/components present in the flow.
geometry *	Integer	Not Supported	Geometry configuration of the patch.
alter_patch(i) *	Logical	Not Supported	Alter the $i$-th patch.
x[y,z]_centroid *	Real	Not Supported	Centroid of the applied geometry in the $[x,y,z]$-direction.
length_x[y,z] *	Real	Not Supported	Length, if applicable, in the $[x,y,z]$-direction.
radius *	Real	Not Supported	Radius, if applicable, of the applied geometry.
smoothen *	Logical	Not Supported	Smoothen the applied patch.
smooth_patch_id *	Integer	Not Supported	A patch with which the applied patch is smoothened.
smooth_coeff *	Real	Not Supported	Smoothen coefficient.
alpha(i) *	Real	Supported	Volume fraction of fluid $i$.
alpha_rho(i) *	Real	Supported	Partial density of fluid $i$.
pres *	Real	Supported	Pressure.
vel(i) *	Real	Supported	Velocity in direction $i$.
tau_e(i) *	Real	Supported	Elastic stresses.
hcid *	Integer	N/A	Hard coded patch id
cf_val *	Real	Supported	Surface tension color function value
model_id	Integer	Not Supported	Index into the stl_models array (geometry 21)

*: These parameters should be prepended with patch_icpp(j)% where $j$ is the patch index.

The Table lists the patch parameters. The parameters define the geometries and physical parameters of fluid components (patch) in the domain at initial condition. Note that the domain must be fully filled with patches. The code outputs error messages when an empty region is left in the domain.

• tau_e(i) is the i-th component of the elastic stress tensor, ordered as tau_xx, tau_xy, tau_yy, tau_xz, tau_yz, and tau_zz. 1D simulation requires tau(1), 2D tau(1:3), and 3D tau(1:6).

Analytical Definition of Primitive Variables

Some parameters, as described above, can be defined by analytical functions in the input file. For example, one can define the following patch:

'patch_icpp(2)%geometry'    : 1,
'patch_icpp(2)%x_centroid'  : 0.25,
'patch_icpp(2)%length_x'    : 9.5,
'patch_icpp(2)%vel(1)'      : 0.,
'patch_icpp(2)%pres'        : 1.,
'patch_icpp(2)%alpha_rho(1)': '1 + 0.1*sin(20*x*pi)',
'patch_icpp(2)%alpha(1)'    : 1.,

where alpha_rho is defined with the 1 + 0.1*sin(20*x*pi) function.

Expressions use Python syntax and are parsed at case load time. Syntax errors and unknown variable or function names are immediate, named errors — they are reported before any Fortran is compiled, with a message identifying the offending expression and listing the available names.

The following variables are available in IC patch expressions:

Shorthand	Expands To	Shorthand	Expands To	Shorthand	Expands To
x	x_cc(i)	lx	The patch's length_x	xc	The patch's x_centroid
y	y_cc(j)	ly	The patch's length_y	yc	The patch's y_centroid
z	z_cc(k)	lz	The patch's length_z	zc	The patch's z_centroid
eps	The patch's epsilon	beta	The patch's beta	radii	The patch's radii
tau_e	The patch's tau_e	r	The patch's radius	pi	\f$\pi\f$ (Fortran constant from m_constants)

where $(i,j,k)$ are the grid-indices of the current cell in each coordinate direction.

The allowed functions are the standard Fortran intrinsics: sin, cos, tan, asin, acos, atan, atan2, sinh, cosh, tanh, exp, log, log10, sqrt, abs, min, max, mod, sign.

Euler's number: bare e is not a variable. Write exp(1.0) or a numeric literal instead (e.g. 2.718281828).

In the example above, the following code is generated:

if (patch_id == 2) then
    q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1 + 0.1 * sin(20 * x_cc(i) * pi)
end if

Hard Coded Patches

Some patch configurations are not adequately handled with the above analytic variable definitions. In this case, a hard coded patch can be used. Hard coded patches can be added by adding additional hard coded patch identifiers to src/common/include/1[2,3]dHardcodedIC.fpp. When using a hard coded patch, the patch_icpp(patch_id)%%hcid must be set to the hard-coded patch id. For example, to add a 2D Hardcoded patch with an id of 200, one would add the following to src/common/include/2dHardcodedIC.fpp

    case(200)
        ! Primitive variables assignment

and use patch_icpp(i)%%hcid = 200 in the input file. Additional variables can be declared in Hardcoded1[2,3]DVariables and used in hardcoded1[2,3]D. As a convention, any hard coded patches that are part of the MFC master branch should be identified as 1[2,3]xx where the first digit indicates the number of dimensions.

The code provides three pre-built patches for dimensional extrusion of initial conditions:

• case(170): Load 1D profile from data files
• case(270): Extrude 1D data to 2D domain
• case(370): Extrude 2D data to 3D domain

Setup: Only requires specifying init_dir and filename pattern via zeros_default. Grid dimensions are automatically detected from the data files. Implementation: All variables and file handling are managed in src/common/include/ExtrusionHardcodedIC.fpp with no manual grid configuration needed. Usage: Ideal for initializing simulations from lower-dimensional solutions, enabling users to add perturbations or modifications to the base extruded fields for flow instability studies.

Parameter Descriptions

• num_patches defines the total number of patches defined in the domain. The number has to be a positive integer.

• num_fluids defines the total number of fluids defined in each of the patches. The number has to be a positive integer.

• patch_icpp(j)%%geometry defines the type of geometry of $j$-th patch by using an integer from 1 to 21. Definition of the patch type for each integer is listed in table Patch Types.

• [x,y,z]_centroid, length_[x,y,z], and/or radius are used to uniquely define the geometry of the patch with given type. Requisite combinations of the parameters for each type are listed in table Patch types.

• patch_icpp(j)%%alter_patch(i) activates alternation of patch(i) with patch(j). For instance, in a 2D simulation, when a cylindrical patch(2) is immersed in a rectangular patch(1):

  • patch_icpp(1)%%geometry = 3
  • patch_icpp(2)%%geometry = 2
  • patch_icpp(2)%%alter_patch(1) = 'T'

• smoothen activates smoothening of the boundary of the patch that alters the existing patch. When smoothening occurs, fluids of the two patches are mixed in the region of the boundary. For instance, in the aforementioned case of the cylindrical patch immersed in the rectangular patch, smoothening occurs when patch_icpp(2)%%smoothen = 'T'. smooth_coeff controls the thickness of the region of smoothening (sharpness of the mixture region). The default value of smooth_coeff is unity. The region of smoothening is thickened with decreasing the value. Optimal choice of the value of smooth_coeff is case-dependent and left to the user.

• patch_icpp(j)%%alpha(i), patch_icpp(j)%%alpha_rho(i), patch_icpp(j)%%pres, and patch_icpp(j)%%vel(i) define for $j$-th patch the void fraction of fluid(i), partial density of fluid(i), the pressure, and the velocity in the $i$-th coordinate direction. These physical parameters must be consistent with fluid material's parameters defined in the next subsection.

• model_id selects the STL/OBJ model for a geometry-21 patch by indexing into the stl_models array. The model file, scaling, and translation, and the inside/outside threshold, are configured on that stl_models entry (see the stl_models section below); a cell is marked inside the model using a winding-number test.

Elliptic Smoothing

Initial conditions in which not all patches support the patch_icpp(j)%%smoothen parameter can still be smoothed by applying iterations of the heat equation to the initial condition. This is enabled by adding 'elliptic_smoothing': "T", and 'elliptic_smoothing_iters': N, to the case dictionary, where N is the number of smoothing iterations to apply.

4. Immersed Boundary Patches {#sec-immersed-boundary-patches}

Parameter	Type	Description
num_ibs	Integer	Number of immersed boundary patches
num_stl_models	Integer	Number of STL/OBJ model entries in the stl_models array
num_particle_clouds	Integer	Number of particle bed specifications to generate immersed boundary patches from
ib_neighborhood_radius	Integer	Parameter that controls the neighborhood size for IB detection.
many_ib_patch_parallelism	Logical	Parallelize over IB patches instead of grid cells (better for many small patches).
geometry	Integer	Geometry configuration of the patch.
x[y,z]_centroid	Real	Centroid of the applied geometry in the [x,y,z]-direction.
length_x[y,z]	Real	Length, if applicable, in the [x,y,z]-direction.
radius	Real	Radius, if applicable, of the applied geometry.
airfoil_id	Integer	Index into ib_airfoil array for NACA airfoil geometry patches.
model_id	Integer	Index into stl_models array for STL/OBJ geometry patches.
slip	Logical	Apply a slip boundary
moving_ibm	Integer	Sets the method used for IB movement.
vel(i)	Real	Initial velocity of the moving IB in the i-th direction.
angular_vel(i)	Real	Initial angular velocity of the moving IB in the i-th direction.
coefficient_of_restitution	Real	A number 0 to 1 describing how elastic IB collisions are
collision_model	Integer	Integer to select the collision model being used for IB collisions.
collision_time	Real	Amount of simulation time used to resolve collisions
ib_coefficient_of_friction	Real	Coefficient of friction used in IB collisions

These parameters should be prepended with patch_ib(j)% where $j$ is the patch index.

STL/OBJ model geometry parameters are set on the stl_models array (indexed by model_id):

Parameter	Type	Description
model_filepath	String	Path to an STL or OBJ file (not all OBJs are supported).
model_scale(i)	Real	Model's scaling factor for component $i$.
model_translate(i)	Real	Model's $i$-th component of translation.
model_threshold	Real	Winding number threshold above which a cell is marked as inside the model.

These parameters should be prepended with stl_models(k)% where $k$ is the model index.

Parameter Descriptions

• geometry defines the type of geometry of an immersed boundary patch with an integer number. Definitions for currently implemented immersed boundary patch types are listed in table Immersed Boundary Patch Type.

• x[y,z]_centroid is the centroid location of the patch in the x[y,z]-direction

• length_x[y,z] is the length of the patch in the x[y,z]-direction.

• radius is the radius to be used for circular patches.

• c, t, p, and m specify the parameters for a NACA airfoil (set on the referenced ib_airfoil entry). m is the maximum camber, p is the location of maximum camber, c is the chord length, and t is the thickness. Additional details on this specification can be found in NACA airfoil.

• slip applies a slip boundary to the surface of the patch if true and a no-slip boundary condition to the surface if false.

• For STL/OBJ geometry (geometry 5 or 12), set model_id to index into the stl_models array and specify model_filepath, model_scale, model_translate, and model_threshold on that entry.

• moving_ibm sets the method by which movement will be applied to the immersed boundary. Using 0 will result in no movement. Using 1 will result 1-way coupling where the boundary moves at a constant rate and applied forces to the fluid based upon its own motion. In 1-way coupling, the fluid does not apply forces back onto the IB. Using 2 will result in 2-way coupling, where the boundary pushes on the fluid and the fluid pushes back on the boundary via pressure and viscous forces. If external forces are applied, the boundary will also experience those forces.

• vel(i) is the initial linear velocity of the IB in the x, y, z direction for i=1, 2, 3. When moving_ibm equals 2, this velocity is just the starting speed of the object, which will then accelerate due to external forces. If moving_ibm equals 1, then this is constant if it is a number, or can be described analytically with an expression.

• angular_vel(i) is the initial angular velocity of the IB about the x, y, z axes for i=1, 2, 3 in radians per second. When moving_ibm equals 2, this rotation rate is just the starting rate of the object, which will then change due to external torques. If moving_ibm equals 1, then this is constant if it is a number, or can be described analytically with an expression.

  Moving-IB analytic expressions use the same Python syntax and error-reporting as IC patch expressions (see the "Analytical Definition of Primitive Variables" section above). Available variables: x (x_cc(i)), y (y_cc(j)), z (z_cc(k)), t (current simulation time), and r (the IB patch radius). The same intrinsic functions and pi constant apply; bare e is not available.

• coefficient_of_restitution is a number from 0 (exclusive) to 1 (inclusive) describing how elastic IB collisions are. 0 is for perfectly inelastic collisions while 1 is for perfectly elastic collisions.

• collision_model is an integer to select the collision model being used for IB collisions. Using 0 disables collisions and collision checking. 1 enables the soft-sphere collision model, where all IBs must be circles or sphere and those IBs can collide with each other as well as walls.

• collision_time is approximately the amount of simulation time used to resolve collisions. This is handled by modifying the spring constant used to apply collision forces.

• ib_coefficient_of_friction is the coefficient of friction used in IB collisions.

• ib_neighborhood_radius controls the size of the neighborhood size. This value defaults to 1, which indicates that any given rank is aware of IBs up to 1 ranks away. This parameter is required to strong-scale a case when IBs eventually grow to be larger than one full processor domain wide.

Particle Clouds

A particle cloud is a compact specification of a bed of identical circular (2D) or spherical (3D) immersed boundaries; each cloud is expanded into individual patch_ib particles at startup. Set num_particle_clouds to the number of beds and prepend the parameters below with particle_cloud(j)% where $j$ is the cloud index.

Parameter	Type	Description
x[y,z]_centroid	Real	Centre of the cloud region in the [x,y,z]-direction.
length_x[y,z]	Real	Extent of the cloud region in the [x,y,z]-direction.
num_particles	Integer	Number of particles to place in the region.
radius	Real	Radius of every particle in the cloud.
mass	Real	Mass of every particle in the cloud.
min_spacing	Real	Minimum surface-to-surface gap between particles (centres are 2*radius + min_spacing apart).
moving_ibm	Integer	Motion flag applied to every particle (see patch_ib(j)%%moving_ibm).
seed	Integer	Random seed for reproducible placement (used by packing_method = 1).
packing_method	Integer	Algorithm used to place the particles.

• packing_method selects how the num_particles are positioned within the cloud region:
  • 1 (rejection sampling) draws random positions and rejects any that violate min_spacing, producing a disordered bed. seed makes the placement reproducible.
  • 2 (lattice) places the particles on the optimally dense lattice for the geometry — a triangular lattice in 2D and a face-centered cubic lattice in 3D. The lattice spacing is derived from the particle density (num_particles over the region area/volume); if that spacing is below the required 2*radius + min_spacing, the region is too dense and the run aborts.

5. Fluid Material's {#sec-fluid-materials}

Parameter	Type	Description
gamma	Real	Stiffened-gas parameter \f$\Gamma\f$ of fluid.
pi_inf	Real	Stiffened-gas parameter \f$\Pi_\infty\f$ of fluid.
Re(1) *	Real	Shear viscosity of fluid.
Re(2) *	Real	Volume viscosity of fluid.
cv **	Real	Sffened-gas parameter $c_v$ of fluid.
qv **	Real	Stiffened-gas parameter $q$ of fluid.
qvp **	Real	Stiffened-gas parameter $q'$ of fluid.
sigma	Real	Surface tension coefficient
G	Real	Shear modulus of solid.

Fluid material's parameters. All parameters except for sigma should be prepended with fluid_pp(i) where $i$ is the fluid index.

*: Parameters that work only with model_eqns = 2.

**: Parameters that work only with model_eqns = 3.

The table lists the fluid material's parameters. The parameters define material's property of compressible fluids that are used in simulation.

• fluid_pp(i)%%gamma and fluid_pp(i)%%pi_inf define \f$\Gamma\f$ and \f$\Pi\f$ as parameters of $i$-th fluid that are used in stiffened gas equation of state.

• fluid_pp(i)%%Re(1) and fluid_pp(i)%%Re(2) define the shear and volume viscosities of $i$-th fluid, respectively.

When these parameters are undefined, fluids are treated as inviscid. Details of implementation of viscosity in MFC can be found in \cite Coralic15.

• fluid_pp(i)%%cv, fluid_pp(i)%%qv, and fluid_pp(i)%%qvp define $c_v$, $q$, and $q'$ as parameters of $i$-th fluid that are used in stiffened gas equation of state.

• fluid_pp(i)%%G is required for hypoelasticity.

Stored-form parameters: The values gamma, pi_inf, and Re(1)/Re(2) are not the raw physical quantities. MFC expects transformed stored forms:

• gamma = \f$1/(\gamma-1)\f$, not \f$\gamma\f$ itself
• pi_inf = \f$\gamma,\pi_\infty / (\gamma - 1)\f$, not \f$\pi_\infty\f$ itself
• Re(1) = \f$1/\mu\f$ (inverse viscosity), not \f$\mu\f$ itself

Setting gamma = 1.4 for air is a common mistake; the correct value is 1.0 / (1.4 - 1.0) = 2.5. See @ref sec-stored-forms and @ref sec-material-values in the Equations reference for the full table.

6. Simulation Algorithm {#sec-simulation-algorithm}

See @ref equations "Equations" for the mathematical models these parameters control.

Parameter	Type	Description
bc_[x,y,z]%%beg[end]	Integer	Beginning [ending] boundary condition in the $[x,y,z]$-direction (negative integer, see table Boundary Conditions)
bc_[x,y,z]%%vb[1,2,3]‡	Real	Velocity in the (x,1), (y, 2), (z,3) direction applied to bc_[x,y,z]%%beg
bc_[x,y,z]%%ve[1,2,3]‡	Real	Velocity in the (x,1), (y, 2), (z,3) direction applied to bc_[x,y,z]%%end
model_eqns	Integer	Multicomponent model: [1] \f$\Gamma/\Pi_\infty\f$; [2] 5-equation; [3] 6-equation; [4] 4-equation
alt_soundspeed *	Logical	Alternate sound speed and \f$K \nabla \cdot u\f$ for 5-equation model
adv_n	Logical	Solving directly for the number density (in the method of classes) and compute void fraction from the number density
mpp_lim	Logical	Mixture physical parameters limits
mixture_err	Logical	Mixture properties correction
time_stepper	Integer	Runge--Kutta order [1-3]
adap_dt	Logical	Strang splitting scheme with adaptive time stepping
recon_type	Integer	Reconstruction Type: [1] WENO; [2] MUSCL
adap_dt_tol	Real	Tolerance for adaptive time stepping in Strang splitting scheme
adap_dt_max_iters	Integer	Max iteration for adaptive time stepping in Strang splitting scheme
weno_order	Integer	WENO order [1,3,5]
weno_eps	Real	WENO perturbation (avoid division by zero)
mapped_weno	Logical	WENO-M (WENO with mapping of nonlinear weights)
wenoz	Logical	WENO-Z
wenoz_q	Real	WENO-Z power parameter q (only for WENO7)
teno	Logical	TENO (Targeted ENO)
teno_CT	Real	TENO threshold for smoothness detection
null_weights	Logical	Null WENO weights at boundaries
mp_weno	Logical	Monotonicity preserving WENO
muscl_order	Integer	MUSCL order [1,2]
muscl_lim	Integer	MUSCL Slope Limiter: [1] minmod; [2] monotonized central; [3] Van Albada; [4] Van Leer; [5] SUPERBEE
muscl_eps	Real	MUSCL limiter slope-product threshold (default: hard-coded thresholds; set to 0 for textbook behavior)
int_comp	Integer	Interface Compression [0] Off [1] THINC [2] MTHINC (default 0)
flux_lim	Integer	Flux limiter for post-process: [1] minmod; [2] MUSCL; [3] OSPRE; [4] SUPERBEE
ic_eps	Real	Interface compression threshold (default: 1e-4)
ic_beta	Real	Interface compression sharpness parameter (default: 1.6)
riemann_solver	Integer	Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*; [4] HLLD (only for MHD)
low_Mach	Integer	Low Mach number correction for HLLC Riemann solver: [0] None; [1] Pressure (\cite Chen22); [2] Velocity (\cite Thornber08)
avg_state	Integer	Averaged state evaluation method: [1] Roe average*; [2] Arithmetic mean
wave_speeds	Integer	Wave-speed estimation: [1] Direct (\cite Batten97); [2] Pressure-velocity* (\cite Toro09)
weno_Re_flux	Logical	Compute velocity gradient using scalar divergence theorem
weno_avg	Logical	Arithmetic mean of left and right, WENO-reconstructed, cell-boundary values
dt	Real	Time step size
t_step_start	Integer	Simulation starting time step
t_step_stop	Integer	Simulation stopping time step
t_step_save	Integer	Frequency to output data
t_step_print	Integer	Frequency to print the current step number to standard output (default 1)
cfl_adap_dt	Logical	CFL based adaptive time-stepping
cfl_const_dt	Logical	CFL based non-adaptive time-stepping
cfl_dt	Logical	Enable CFL-based time stepping
cfl_target	Real	Specified CFL value
n_start	Integer	Save file from which to start simulation
t_save	Real	Time duration between data output
t_stop	Real	Simulation stop time
surface_tension	Logical	Activate surface tension
viscous	Logical	Activate viscosity
hypoelasticity	Logical	Activate hypoelasticity*
pre_stress	Logical	Enable pre-stress initialization for hypoelasticity
igr	Logical	Enable solution via information geometric regularization (IGR) \cite Cao24
igr_order	Integer	Order of reconstruction for IGR [3,5]
alf_factor	Real	Alpha factor for IGR entropic pressure (default 10)
igr_pres_lim	Logical	Limit IGR pressure to avoid negative values (default F)
igr_iter_solver	Integer	Solution method for IGR elliptic solve [1] Jacobi [2] Gauss-Seidel
num_igr_iters	Integer	Number of iterations for for the IGR elliptic solve (default 2)
num_igr_warm_start_iters	Integer	Number of iterations for the IGR elliptic solve at the first time step (default 50)

• * Options that work only with model_eqns = 2.
• † Options that work only with cyl_coord = 'F'.
• ‡ Options that work only with bc_[x,y,z]%[beg,end] = -15 and/or bc_[x,y,z]%[beg,end] = -16.

The table lists simulation algorithm parameters. The parameters are used to specify options in algorithms that are used to integrate the governing equations of the multi-component flow based on the initial condition. Models and assumptions that are used to formulate and discretize the governing equations are described in \cite Wilfong26 and \cite Bryngelson21. Details of the simulation algorithms and implementation of the WENO scheme can be found in \cite Coralic15.

• bc_[x,y,z]%[beg,end] specifies the boundary conditions at the beginning and the end of domain boundaries in each coordinate direction by a negative integer from -1 through -16. See table Boundary Conditions for details. Boundary condition patches can be used with non-characteristic boundary conditions. Their use is detailed in Boundary Condition Patches.

• bc_[x,y,z]%%vb[1,2,3] specifies the velocity in the (x,1), (y,2), (z,3) direction applied to bc_[x,y,z]%%beg when using bc_[x,y,z]%%beg = -16. Tangential velocities require viscosity, weno_avg = T, and bc_[x,y,z]%%beg = -16 to work properly. Normal velocities require bc_[x,y,z]%%end = -15 or bc_[x,y,z]%%end = -16 to work properly.

• bc_[x,y,z]%%ve[1,2,3] specifies the velocity in the (x,1), (y,2), (z,3) direction applied to bc_[x,y,z]%%beg when using bc_[x,y,z]%%end = -16. Tangential velocities require viscosity, weno_avg = T, and bc_[x,y,z]%%end = 16 to work properly. Normal velocities require bc_[x,y,z]%%end = -15 or bc_[x,y,z]%%end = -16 to work properly.

• model_eqns specifies the choice of the multi-component model that is used to formulate the dynamics of the flow using integers from 1 through 3. model_eqns = 1, 2, and 3 correspond to \f$\Gamma\f$-\f$\Pi_\infty\f$ model (\cite Johnsen08), 5-equation model (\cite Allaire02), and 6-equation model (\cite Saurel09), respectively. The difference of the two models is assessed by (\cite Schmidmayer20). Note that some code parameters are only compatible with 5-equation model.

• alt_soundspeed activates the source term in the advection equations for the volume fractions, \f$K\nabla\cdot \underline{u}\f$, that regularizes the speed of sound in the mixture region when the 5-equation model is used. The effect and use of the source term are assessed by \cite Schmidmayer20.

• adv_n activates the direct computation of number density by the Riemann solver instead of computing number density from the void fraction in the method of classes.

• mpp_lim activates correction of solutions to avoid a negative void fraction of each component in each grid cell, such that \f$\alpha_i>\varepsilon\f$ is satisfied at each time step.

• mixture_err activates correction of solutions to avoid imaginary speed of sound at each grid cell.

• time_stepper specifies the order of the Runge-Kutta (RK) time integration scheme that is used for temporal integration in simulation, from the 1st to 5th order by corresponding integer. Note that time_stepper = 3 specifies the total variation diminishing (TVD), third order RK scheme (\cite Gottlieb98).

• adap_dt activates the Strang operator splitting scheme which splits flux and source terms in time marching, and an adaptive time stepping strategy is implemented for the source term. It requires bubbles_euler = 'T', polytropic = 'T', adv_n = 'T' and time_stepper = 3. Additionally, it can be used with bubbles_lagrange = 'T' and time_stepper = 3. adap_dt_tol and adap_dt_max_iters are 1e-4 and 100, respectively, by default.

• weno_order specifies the order of WENO scheme that is used for spatial reconstruction of variables by an integer of 1, 3, 5, and 7, that correspond to the 1st, 3rd, 5th, and 7th order, respectively.

• weno_eps specifies the lower bound of the WENO nonlinear weights. It is recommended to set weno_eps to $10^{-6}$ for WENO-JS, and to $10^{-40}$ for other WENO variants.

• mapped_weno activates the WENO-M scheme in place of the default WENO-JS scheme (\cite Henrick05). WENO-M a variant of the WENO scheme that remaps the nonlinear WENO-JS weights by assigning larger weights to non-smooth stencils, reducing dissipation compared to the default WENO-JS scheme, at the expense of higher computational cost. Only one of mapped_weno, wenoz, and teno can be activated.

• wenoz activates the WENO-Z scheme in place of the default WENO-JS scheme (\cite Borges08). WENO-Z is a variant of the WENO scheme that further reduces the dissipation compared to the WENO-M scheme. It has similar computational cost to the WENO-JS scheme.

• wenoz_q specifies the power parameter q used in the WENO-Z scheme. It controls how aggressively the smoothness coefficients scale the weights. A higher value of wenoz_q increases the sensitivity to smoothness, improving stability but worsening numerical dissipation. For WENO3 and WENO5, q=1 is fixed, so wenoz_q must not be set. For WENO7, wenoz_q can be set to 2, 3, or 4.

• teno activates the TENO scheme in place of the default WENO-JS scheme (\cite Fu16). TENO is a variant of the ENO scheme that is the least dissipative, but could be less robust for extreme cases. It uses a threshold to identify smooth and non-smooth stencils, and applies optimal weights to the smooth stencils. Only available for weno_order = 5 and 7. Requires teno_CT to be set. Does not support grid stretching.

• teno_CT specifies the threshold for the TENO scheme. This dimensionless constant, also known as $C_T$, sets a threshold to identify smooth and non-smooth stencils. Larger values make the scheme more robust but also more dissipative. A recommended value for teno_CT is 1e-6. When adjusting this parameter, it is recommended to try values like 1e-5 or 1e-7 for TENO5. A smaller value can be used for TENO7.

• null_weights activates nullification of the nonlinear WENO weights at the buffer regions outside the domain boundaries when the Riemann extrapolation boundary condition is specified (bc_[x,y,z]%%beg[end]} = -4).

• mp_weno activates monotonicity preservation in the WENO reconstruction (MPWENO) such that the values of reconstructed variables do not reside outside the range spanned by WENO stencil (\cite Balsara00; \cite Suresh97).

• muscl_order specifies the order of the MUSCL scheme that is used for spatial reconstruction of variables by an integer of 1, or 2, that corresponds to the 1st, and 2nd order respectively. When using muscl_order = 2, muscl_lim must be defined.

• muscl_lim specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5. muscl_lim = 1, 2, 3, 4, and 5 correspond to minmod, monotonized central, Van Albada, Van Leer, and SUPERBEE, respectively.

• muscl_eps controls the slope-product activation threshold for all MUSCL limiters. When not set (default), the threshold is 1e-9 for minmod/MC, and 1e-6 for others. Setting muscl_eps = 0 gives textbook limiter behavior where limiters activate whenever both slopes have the same sign.

• int_comp activates interface compression using [1] THINC or [2] MTHINC (default off) used in variable reconstruction, with control parameters (ic_eps, and ic_beta).

• riemann_solver specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 4. riemann_solver = 1, 2, and 3 correspond to HLL, HLLC, and Exact Riemann solver, respectively (\cite Toro09). riemann_solver = 4 is only for MHD simulations. It resolves 5 of the full seven-wave structure of the MHD equations (\cite Miyoshi05).

• low_Mach specifies the choice of the low Mach number correction scheme for the HLLC Riemann solver. low_Mach = 0 is default value and does not apply any correction scheme. low_Mach = 1 and 2 apply the anti-dissipation pressure correction method (\cite Chen22) and the improved velocity reconstruction method (\cite Thornber08). This feature requires model_eqns = 2 or 3. low_Mach = 1 works for riemann_solver = 1 and 2, but low_Mach = 2 only works for riemann_solver = 2.

• avg_state specifies the choice of the method to compute averaged variables at the cell-boundaries from the left and the right states in the Riemann solver by an integer of 1 or 2. avg_state = 1 and 2 correspond to Roe- and arithmetic averages, respectively.

• wave_speeds specifies the choice of the method to compute the left, right, and middle wave speeds in the Riemann solver by an integer of 1 and 2. wave_speeds = 1 and 2 correspond to the direct method (\cite Batten97), and indirect method that approximates the pressures and velocity (\cite Toro09), respectively.

• weno_Re_flux activates the scalar divergence theorem in computing the velocity gradients using WENO-reconstructed cell boundary values. If this option is false, velocity gradient is computed using finite difference scheme of order 2 which is independent of the WENO order.

• weno_avg it activates the arithmetic average of the left and right, WENO-reconstructed, cell-boundary values. This option requires weno_Re_flux to be true because cell boundary values are only utilized when employing the scalar divergence method in the computation of velocity gradients.

• surface_tension activates surface tension when set to 'T'. Requires sigma to be set and num_fluids. The color function in each patch should be assigned such that patch_icpp(i)%%cf_val = 1 in patches where patch_icpp(i)%%alpha = 1 - eps and patch_icpp(i)%%cf_val = 0 in patches where patch_icpp(i)%%alpha = eps.

• viscous activates viscosity when set to 'T'. Requires Re(1) and Re(2) to be set.

• hypoelasticity activates elastic stress calculations for fluid-solid interactions. Requires G to be set in the fluid material's parameters.

Boundary Condition Patches {#boundary-condition-patches}

Parameter	Type	Description
num_bc_patches	Integer	Number of boundary condition patches
dir*	Integer	Direction of the boundary patch. [1]: x; [2]: y; [3]: z
loc*	Integer	Location of the patch in the domain
type*	Integer	The geometry of the patch. [1]: Line [2]: Circle [3]: Rectangle
x[y,z]_centroid*	Real	Centroid of the boundary patch in the x[y,z]-direction
length_x[y,z]*	Real	Length of the boundary patch in the x[y,z]-direction
radius*	Real	Radius of the boundary patch
*: These parameters should be prepended with patch_bc(j)% where $j$ is the patch index.

Boundary condition patches can be used with the following boundary condition types:

• -2 reflective
• -3 ghost cell extrapolation
• -15 slip wall
• -16 no-slip wall
• -17 Dirichlet

Line segments along each domain edge are supported for 2D simulations. Squares and circles on each face are supported for 3D simulations.

Constant Time-Stepping

• dt specifies the constant time step size used in the simulation. The value of dt needs to be sufficiently small to satisfy the Courant-Friedrichs-Lewy (CFL) condition.

• t_step_start and t_step_stop define the time steps at which the simulation starts and ends.

t_step_save is the time step interval for data output during simulation. To newly start the simulation, set t_step_start = 0. To restart the simulation from $k$-th time step, set t_step_start = k; see @ref running "Restarting Cases".

Adaptive Time-Stepping

• cfl_adap_dt enables adaptive time stepping with a constant CFL when true

• cfl_const_dt enables constant dt time-stepping where dt results in a specified CFL for the initial condition

• cfl_target specifies the target CFL value

• n_start specifies the save file to start at

• t_save specifies the time interval between data output during the simulation

• t_stop specifies at what time the simulation should stop

To newly start the simulation, set n_start = 0. To restart the simulation from $k$-th time step, see @ref running "Restarting Cases".

7. Formatted Output {#sec-formatted-output}

Parameter	Type	Description
format	Integer	Output format. [1]: Silo-HDF5; [2] Binary
precision	Integer	[1] Single; [2] Double
parallel_io	Logical	Parallel I/O
file_per_process	Logical	Whether or not to write one IO file per process
cons_vars_wrt	Logical	Write conservative variables
prim_vars_wrt	Logical	Write primitive variables
load_weight_wrt	Logical	Write per-cell load-weight diagnostic field
sfc_partition_wrt	Logical	Report SFC-weighted load-balance partition
rank_time_wrt	Logical	Report per-rank RHS compute-time imbalance (max/mean)
load_balance	Logical	(Experimental/diagnostic) Weighted static Cartesian decomposition at init (requires parallel_io = T, >1 rank). Measured gain is small on CPU (~5%) and can be slower on GPU due to the occupancy floor; equal decomposition is near-optimal for uniform-cost workloads.
amr	Logical	(Experimental) Enable block-structured AMR: a 2:1 refined level-1 block with gradient-based dynamic regrid, optional dt/2 subcycling, and conservative coupling with refluxing. Requires WENO reconstruction, SSP-RK3, model_eqns=2 or 3; num_fluids > 1 requires mpp_lim; supports physical viscosity.
amr_block_beg(i)	Integer	Refined-block start cell index in direction $i$ (level-0 index space)
amr_block_end(i)	Integer	Refined-block end cell index in direction $i$ (level-0 index space)
amr_regrid_int	Integer	Steps between AMR regrid events (0 = static block)
amr_tag_eps	Real	Relative density-gradient threshold for AMR refinement tagging (default 0.1)
amr_buf	Integer	Coarse-cell padding around tagged cells when regridding (default 3)
amr_subcycle	Logical	Advance the coarse level at the case dt and the fine level at dt/2 (two substeps; Berger-Colella refluxing). Requires amr; incompatible with cfl_dt.
amr_max_blocks	Integer	Number of fixed refined-block slots preallocated (each max-block sized; ~N x device memory); must be >= 1 (default 4)
amr_cluster_eff	Real	Berger-Rigoutsos min tag efficiency a clustered block box reaches before splitting stops; must satisfy 0 < eff <= 1 (default 0.7)
hybrid_weno	Logical	Use linear-optimal reconstruction in smooth cells, full WENO only at flagged discontinuities (requires WENO reconstruction)
hybrid_weno_eps	Real	Smoothness threshold for hybrid WENO shock flagging; must be > 0 (default 1e-2)
hybrid_riemann	Logical	Use a cheap central/Rusanov flux in smooth cells, full HLLC only at flagged discontinuities (requires HLLC, 5eq/6eq)
hybrid_smooth_flux	Integer	Smooth-region flux for hybrid Riemann: 1 = central, 2 = Rusanov (default 2)
partition_tile_size	Integer	Tile side for the SFC partitioner (default 8)
alpha_rho_wrt(i)	Logical	Add the partial density of the fluid $i$ to the database |
rho_wrt	Logical	Add the mixture density to the database
mom_wrt(i)	Logical	Add the $i$-direction momentum to the database
vel_wrt(i)	Logical	Add the $i$-direction velocity to the database
E_wrt	Logical	Add the total energy to the database
pres_wrt	Logical	Add the pressure to the database
alpha_wrt(i)	Logical	Add the volume fraction of fluid $i$ to the database
gamma_wrt	Logical	Add the specific heat ratio function to the database
heat_ratio_wrt	Logical	Add the specific heat ratio to the database
ib_state_wrt	Logical	Parameter to handle writing IB state on saves and outputting the state as a point mesh to SILO files.
pi_inf_wrt	Logical	Add the liquid stiffness function to the database
pres_inf_wrt	Logical	Add the liquid stiffness to the formatted database
c_wrt	Logical	Add the sound speed to the database
omega_wrt(i)	Logical	Add the $i$-direction vorticity to the database
schlieren_wrt	Logical	Add the numerical schlieren to the database
qm_wrt	Logical	Add the Q-criterion to the database
liutex_wrt	Logical	Add the Liutex to the database
cf_wrt	Logical	Write color function field
chem_wrt_T	Logical	Write temperature field for chemistry output
fft_wrt	Logical	Enable FFT output
sim_data	Logical	Write interface and energy data files (post_process)
integral_wrt	Logical	Write integral data
num_integrals	Integer	Number of integral regions
down_sample	Logical	Enable output downsampling
fd_order	Integer	Order of finite differences for computing the vorticity and the numerical Schlieren function [1,2,4]
schlieren_alpha(i)	Real	Intensity of the numerical Schlieren computed via alpha(i)
probe_wrt	Logical	Write the flow chosen probes data files for each time step
num_probes	Integer	Number of probes
probe(i)%[x,y,z]	Real	Coordinates of probe $i$
output_partial_domain	Logical	Output part of the domain
[x,y,z]_output%%beg	Real	Beginning of the output domain in the [x,y,z]-direction
[x,y,z]_output%%end	Real	End of the output domain in the [x,y,z]-direction
lag_txt_wrt	Logical	Write Lagrangian bubble data to .dat files
lag_header	Logical	Write header to Lagrangian bubble .dat files
lag_db_wrt	Logical	Write Lagrangian bubble data to silo/hdf5 database files
lag_id_wrt	Logical	Add the global bubble idea to the database file
lag_pos_wrt	Logical	Add the bubble position to the database file
lag_pos_prev_wrt	Logical	Add the previous bubble position to the database file
lag_vel_wrt	Logical	Add the bubble translational velocity to the database file
lag_rad_wrt	Logical	Add the bubble radius to the database file
lag_rvel_wrt	Logical	Add the bubble radial velocity to the database file
lag_r0_wrt	Logical	Add the bubble initial radius to the database file
lag_rmax_wrt	Logical	Add the bubble maximum radius to the database file
lag_rmin_wrt	Logical	Add the bubble minimum radius to the database file
lag_dphidt_wrt	Logical	Add the bubble subgrid velocity potential to the database file
lag_pres_wrt	Logical	Add the bubble pressure to the database file
lag_mv_wrt	Logical	Add the bubble vapor mass to the database file
lag_mg_wrt	Logical	Add the bubble gas mass to the database file
lag_betaT_wrt	Logical	Add the bubble heat flux model coefficient to the database file
lag_betaC_wrt	Logical	Add the bubble mass flux model coefficient to the database file

The table lists formatted database output parameters. The parameters define variables that are outputted from simulation and file types and formats of data as well as options for post-processing.

• format specifies the choice of the file format of data file outputted by MFC by an integer of 1 and 2. format = 1 and 2 correspond to Silo-HDF5 format and binary format, respectively. Both formats are supported by ./mfc.sh viz (see @ref visualization "Flow Visualization"). Silo-HDF5 requires the h5py Python package; binary has no extra dependencies.

• precision specifies the choice of the floating-point format of the data file outputted by MFC by an integer of 1 and 2. precision = 1 and 2 correspond to single-precision and double-precision formats, respectively.

• parallel_io activates parallel input/output (I/O) of data files. It is highly recommended to activate this option in a parallel environment. With parallel I/O, MFC inputs and outputs a single file throughout pre-process, simulation, and post-process, regardless of the number of processors used. Parallel I/O enables the use of different numbers of processors in each process (e.g., simulation data generated using 1000 processors can be post-processed using a single processor).

• file_per_process deactivates shared file MPI-IO and activates file per process MPI-IO. The default behavior is to use a shared file. File per process is useful when running on >10K ranks. If file_per_process is true, then pre_process, simulation, and post_process must be run with the same number of ranks.

• cons_vars_wrt and prim_vars_wrt activate the output of conservative and primitive state variables into the database.

• [variable's name]_wrt activates the output of each specified variable into the database.

• schlieren_alpha(i) specifies the intensity of the numerical Schlieren of $i$-th component. It must be specified for every fluid when schlieren_wrt is enabled.

• fd_order specifies the order of the finite difference scheme used to compute the vorticity from the velocity field and the numerical schlieren from the density field using an integer of 1, 2, and 4. fd_order = 1, 2, and 4 correspond to the first, second, and fourth-order finite difference schemes.

• probe_wrt activates the output of state variables at coordinates specified by probe(i)%[x;y,z].

• ib_state_wrt is used to trigger post-processing of the IB state to be written out as a point mesh in the SILO files. When no IBs are moving, it also triggers force and torque calculation so that those values may be written to the output state files.

• output_partial_domain activates the output of part of the domain specified by [x,y,z]_output%%beg and [x,y,z]_output%%end. This is useful for large domains where only a portion of the domain is of interest. It is not supported when precision = 1 and format = 1. It also cannot be enabled with flux_wrt, heat_ratio_wrt, pres_inf_wrt, c_wrt, omega_wrt, ib, schlieren_wrt, qm_wrt, or 'liutex_wrt'.

7.1. Adaptive Mesh Refinement (AMR) {#sec-amr}

MFC supports block-structured AMR (Experimental) via up to amr_max_blocks 2:1 refined level-1 blocks that coexists with the base-level solve. The fine block is initialized from the base grid by piecewise-linear interpolation and remains continuously coupled to the base solve through conservative ghost-cell exchange and flux refluxing at the coarse–fine interface.

Restrictions. AMR requires WENO reconstruction (recon_type = 1, any order), SSP-RK3 time-stepping (time_stepper = 3), and the 5- or 6-equation model (model_eqns = 2 or 3; for 6-eq the per-stage pressure relaxation also runs on each fine block). Multiple fluids (num_fluids > 1) are supported and additionally require mpp_lim, whose volume-fraction clamp+renormalize maintains coarse/fine alpha consistency; the per-fluid masses are refluxed exactly, and volume fractions are prolonged with a sum-preserving closure (fine-level volume fractions sum to one by construction). Physical viscosity (viscous = T) is supported: the viscous stress/work travels through the momentum- and energy-equation source fluxes, which are captured into the same coarse–fine flux registers as the advective fluxes, so the interface is refluxed against the matched total (advective + viscous) flux and energy — including viscous work — is conserved. Fine-ghost velocity gradients at the coarse–fine boundary are taken from the conservative-linear prolongation of the coarse state (no special gradient reconstruction); that interface inconsistency is bounded and conservation is enforced by the flux-register matching. The density-gradient regrid tagger does not sense shear or boundary layers well, so viscous features may need a static or generously buffered block (error-estimator taggers are future work). Euler-Euler bubbles (bubbles_euler = T) are supported, including non-polytropic (polytropic = F) and polydisperse (nb > 1) configurations: the bubble moments — radius, velocity, and, for non-polytropic, partial pressure and vapor mass, per R0 bin — are all flux-based conserved variables refluxed through the same registers, so no separate side-state is carried on the fine level. Prolongation floors every positive moment (radius, and the non-polytropic partial-pressure / vapor-mass moments) while leaving the signed velocity moment free, so the reconstructed radius, number density, internal pressure, and vapor mass stay non-negative (realizability). QBMM (qbmm = T) is supported for the polytropic model: each R0 bin's bivariate six-moment set lives entirely in the conserved variables (the pb/mv quadrature arrays are inert stubs when polytropic = T), and the whole moment block is prolonged piecewise-constant so every fine/ghost cell inherits the coarse cell's realizable moment set (the CHyQMOM inversion needs the radius variance c20 = m20/m00 - (m10/m00)^2 to stay positive, which a per-component minmod slope could break); the moments still reflux and restrict on the standard conservative path. Non-polytropic QBMM (polytropic = F) is fully supported: each block carries its own per-quadrature-node internal pressure and vapor mass (pb/mv), prolonged piecewise-constant for realizability, advanced with the block's own rhs scratch, and restricted back with the moments; dynamic regrid and amr_subcycle are both supported. Phase change (relax) is supported: the cell-local, mass/energy-conserving relaxation runs on the fine solution before restriction (matching the coarse once-per-step timing). Chemistry (chemistry = T) is supported for reactions and advection: the species partial densities are flux-based conserved variables refluxed through the same registers, the cell-local reaction source runs on the fine block through the shared RHS (matching the coarse per-stage timing), and prolongation rescales the species so their sum equals the continuity density (sum(Y_k) = 1, Y_k >= 0 on the fine level by construction). Chemistry AMR runs single- and multi-rank: the fine block's cons->prim conversion widens over the ghost shell, so the temperature (the reacting-EOS Newton guess) is halo-exchanged with the coarse state at rank seams (mirroring the diffusion path) — without it the seam-ghost guess is uninitialized and the conversion diverges to NaN. Species mass diffusion (chem_params%%diffusion = T) is also supported: the mixture-averaged species mass fluxes (and the thermal-conduction + enthalpy energy flux) travel through the source-flux array and are captured into the same coarse–fine registers as the advective fluxes — like the viscous stress fluxes — so element mass and total energy conserve across the block boundary through refluxed, subcycled, and regridded advances. Static immersed boundaries (ib = T) are supported: each fine block carries its own fine-grid IB state (markers, ghost points, levelset, image points, interpolation coefficients) computed from the body geometry at fine resolution once at initialization, and the fine advance applies the ghost-cell IB state correction on the block after each RK update (mirroring the coarse per-stage timing). A fixed body placed inside a static block is thus resolved on the refined level. The IB forcing is non-conservative by construction (the ghost-cell method injects mass/momentum/energy at the body), so the conservation defect is nonzero in the body region while the flux reflux still conserves to machine precision away from it. A body in prescribed motion (moving_ibm = 1) is also supported: the fine block's IB markers/ghost points are rebuilt each fine RK substage at the body's sub-time position (the same linear time interpolation the subcycle applies to the fluid ghosts), so the refined body tracks its prescribed trajectory. Supports one or more non-STL bodies, static or in prescribed motion; with dynamic regrid every candidate box expands to fully contain each body at its live position plus a margin, the fine IB state is rebuilt from geometry after each regrid, and a per-substage guard aborts if a moving body reaches its block boundary between regrids (reduce amr_regrid_int or increase amr_buf). Force-driven motion (moving_ibm = 2) and STL geometry are gated pending validation. Under MPI a body contained within one rank's subdomain is bit-exact across decompositions; a body spanning a rank seam is rejected at startup (the fine-IB image-point stencil across the seam is not yet decomposition-exact), so keep the body inside a single rank's subdomain (use fewer ranks or reposition it). Lagrangian bubbles are supported with the cloud excluded from fine blocks: two-way coupling lives on the coarse grid, regrid suppresses tags and clips boxes around the cloud's padded bounding box, EL volume fractions prolong without the sum-to-one closure (their sum is the local liquid fraction, not 1), and a per-stage guard aborts if the cloud reaches an active block (reduce amr_regrid_int or increase amr_buf). The IGR solver is supported with restriction-only coarse/fine coupling: the fine block runs its own fixed-iteration sigma solve seeded and Dirichlet-bounded by the converged coarse sigma; seam conservation is truncation-order (no reflux capture from the fused IGR flux kernels), free-stream preservation is exact, and amr_subcycle is gated under IGR. AMR is incompatible with surface tension, 3D cylindrical coordinates (2D axisymmetric IS supported), 2D/3D MHD (measured: the coarse/fine seam is a continuous div(B) source that GLM cleaning cannot remove; 1D MHD/RMHD IS supported since div(B) = 0 by construction there), hyperelasticity, and Riemann-extrapolation boundaries (bc = -4). active_box is supported (single-rank): blocks must sit strictly inside the growing active window (init abort + regrid clamp), and the fine advance treats its whole block as active. hybrid_weno/hybrid_riemann are supported: each level recomputes the smoothness sensor over its own (swapped) bounds every RHS call. Nonuniform grids ARE supported (grid stretching and the axisymmetric axis half-cell): the fine ghost-shell coordinates extend by exact parent-cell bisection and the spacing-dependent WENO coefficients are recomputed for the active grid on every block swap/restore, armed automatically when the grid is detected nonuniform at startup. Acoustic sources are supported on the coarse grid only: the support must not overlap the initial block (startup abort) and dynamic regrid keeps its boxes clear of the support. Multi-rank runs are supported: the fine level mirrors the base decomposition (each rank holds the fine cells covering the block's intersection with its own subdomain), so the block may span rank boundaries and move freely across them under dynamic regrid. The block may cover at most about half of any rank's subdomain per dimension (the fine advance reuses the rank-local solver scratch).

AMR + surface tension (unsupported). Surface tension (surface_tension = T) is prohibited under AMR. What works: the capillary contribution is a face flux captured into the same coarse–fine registers as the advective flux, so it is refluxed conservatively — conservation is structural (mass and energy defects stay at machine precision regardless of the fine-side treatment). What fails: the capillary stress is normalized (∝ 1/|∇c|), so it depends on the interface-normal direction, not the gradient magnitude. Across a 2:1 coarse/fine boundary the conservative-linearly-prolonged fine ghost color cannot reproduce the coarse solver's interface normal, producing a growing spurious seam current. Every fine-block fix attempted failed: opening the capillary reflux gate alone (~540x baseline seam velocity, exponential), a smoothstep ramp suppressing the fine capillary force near the seam (~27x, bounded-linear, width-invariant), and a coarse-spacing gradient blend of the prolonged color (~556x, growing) — all leave a force imbalance from the inconsistent interface normal rather than a curvature spike that can be damped. What might fix it: capturing the native coarse-computed capillary force Ω in a per-block band during the coarse RHS, prolonging it to the fine boundary layer, and blending the force there — large and uncertain, and the diffuse-interface 2:1 normal inconsistency may be fundamental.

Static vs. dynamic block. Setting amr_regrid_int = 0 fixes the block at the initial amr_block_beg/amr_block_end position for the entire run (useful for convergence studies or GPU correctness testing). Setting amr_regrid_int > 0 triggers dynamic regrid every that many coarse steps: cells whose normalized density gradient exceeds amr_tag_eps are tagged, then clustered by a Berger–Rigoutsos recursive bisection into a list of separated block boxes (each grown by amr_buf coarse cells of buffer padding). Boxes whose padded extents would come within a ghost-cell buffer width of each other are merged, so separated features get their own refined box while nearby ones stay a single box (guaranteeing no fine–fine adjacency). Splitting stops once a box's tag efficiency (tagged/total cells) reaches amr_cluster_eff; the number of blocks is capped at amr_max_blocks. A positive amr_tag_eps and amr_buf >= 1 are required whenever regridding is active.

Subcycling. amr_subcycle = T enables Berger–Colella dt/2 subcycling: the coarse level advances one full step at the case dt, while the fine level takes two half-steps at dt/2 with time-interpolated ghost values at the intermediate stage. Accumulated fine-level fluxes are applied back to the coarse level (reflux correction) after each coarse step. amr_subcycle is incompatible with cfl_dt (variable time step) and requires amr = T.

Block slots. amr_max_blocks (default 4) sets the number of fixed refined-block slots preallocated for the run. Each slot is sized to the maximum block extent, so N slots require roughly N times the device memory of a single block; the goal is the compute win of refining separated features independently, and memory efficiency (compact per-block pools) is a follow-up. Dynamic regrid clusters the tagged cells into up to amr_max_blocks separated boxes (amr_cluster_eff sets the min tag efficiency each box reaches before splitting stops).

Restart. Each save step writes a fine-level AMR restart file alongside the level-0 restart data (whose format is unchanged): the current — possibly regridded — block box and the fine solution, per rank (an amr_fine.dat in each rank's step directory, or a single shared amr_*.dat next to the level-0 MPI-IO restart file when parallel_io is on). Restarting (t_step_start > 0) restores the saved box and fine state seamlessly; it requires the same rank count (and decomposition) as the run that wrote the file, and the same physics configuration — the number of conserved variables, which depends on num_fluids, model_eqns, and the enabled bubble/chemistry models — and aborts with a clear message otherwise. On restart the AMR block geometry (block count and boxes) is read from the AMR restart file, not from the amr_block_beg/amr_block_end case parameters — so editing those parameters for a restart run has no effect. To re-derive the blocks from parameters, start fresh (t_step_start = 0). If the AMR file is absent (e.g., data from an older run), the run proceeds with a warning and re-initializes the fine level by prolongation from the coarse restart data, losing the accumulated fine-level accuracy. Note that level-0 output already contains the restricted (coarse-resolution) fine solution over the block, so existing visualization works unchanged; fine-resolution visualization output is future work.

Parameter	Type	Description
amr	Logical	Enable AMR (see prose above for requirements and restrictions)
amr_block_beg(i)	Integer	Initial refined-block start cell index in direction $i$ (level-0 index space)
amr_block_end(i)	Integer	Initial refined-block end cell index in direction \f$i\f$ (level-0 index space); must satisfy \f$2,(e_i - b_i + 1) - 1 \le N_i\f$
amr_regrid_int	Integer	Coarse steps between regrid events (0 = static block)
amr_tag_eps	Real	Normalized density-gradient threshold for refinement tagging; must be > 0 when amr_regrid_int > 0 (default 0.1)
amr_buf	Integer	Coarse-cell padding around tagged cells; must be >= 1 when amr_regrid_int > 0 (default 3)
amr_subcycle	Logical	Advance fine level at dt/2 (two substeps per coarse step) with Berger–Colella refluxing
amr_max_blocks	Integer	Number of fixed refined-block slots preallocated (each max-block sized; ~N x device memory); must be >= 1 (default 4)
amr_cluster_eff	Real	Berger-Rigoutsos min tag efficiency a clustered block box reaches before splitting stops; must satisfy 0 < eff <= 1 (default 0.7)

8. Acoustic Source {#sec-acoustic-source}

Parameter	Type	Description
acoustic_source	Logical	Acoustic source module activation
num_source	Integer	Number of acoustic sources
acoustic(i)%%support	Integer	Geometry of spatial support for the acoustic source
acoustic(i)%%dipole	Logical	Dipole source activation (optional; default = false -> monopole)
acoustic(i)%%loc(j)	Real	$j$-th coordinate of the point that defines the acoustic source location
acoustic(i)%%pulse	Integer	Acoustic wave form: [1] Sine [2] Gaussian [3] Square [4] Broadband
acoustic(i)%%npulse	Real	Number of pulse cycles
acoustic(i)%%mag	Real	Pulse magnitude
acoustic(i)%%frequency	Real	Sine/Square - Frequency of the acoustic wave (exclusive)
acoustic(i)%%wavelength	Real	Sine/Square - Wavelength of the acoustic wave (exclusive)
acoustic(i)%%gauss_sigma_time	Real	Gaussian - Gaussian pulse time width in terms of sigma (exclusive)
acoustic(i)%%gauss_sigma_dist	Real	Gaussian - Gaussian pulse spatial width in terms of sigma (exclusive)
acoustic(i)%%delay	Real	Time delay of the acoustic wave (optional for %%pulse = 1 or 3; default = 0)
acoustic(i)%%dir	Real	Planer - Direction of acoustic propagation
acoustic(i)%%length	Real	2D/3D Planer - Spatial pulse length
acoustic(i)%%height	Real	3D Planer - Spatial pulse height
acoustic(i)%%foc_length	Real	Transducer - Focal length of the transducer
acoustic(i)%%aperture	Real	Transducer - Aperture of the transducer
acoustic(i)%%num_elements	Integer	Transducer array - Number of transducer elements in a transducer array
acoustic(i)%%element_on	Integer	Transducer array - Element number that is on (optional; default = 0 -> all elements)
acoustic(i)%%element_spacing_angle	Real	2D Transducer array - Spacing angle (in rad) between adjacent transducer elements
acoustic(i)%%element_polygon_ratio	Real	3D Transducer array - Ratio of polygon side length to transducer element radius
acoustic(i)%%rotate_angle	Real	3D Transducer array - Rotation angle of the transducer array (optional; default = 0)
acoustic(i)%%bb_num_freq	integer	Number of frequencies in broadband wave
acoustic(i)%%bb_bandwidth	Real	The bandwidth of each frequency in the broadband wave
acoustic(i)%%bb_lowest_freq	Real	The lower frequency bound of the broadband wave

Details of the transducer acoustic source model can be found in \cite Maeda17.

• acoustic_source activates the acoustic source module.

• num_source defines the total number of source planes by an integer.

• %%support specifies the choice of the geometry of acoustic source distribution. See table Acoustic Supports for details.

• %%dipole changes the default monopole (one-sided) source to a dipole source. It is only available for planar waves.

• %%loc(j) specifies the location of the acoustic source in the $j$-th coordinate direction. For planer support, the location defines midpoint of the source plane. For transducer arrays, the location defines the center of the transducer or transducer array (not the focal point; for 3D it's the tip of the spherical cap, for 2D it's the tip of the arc).

• %%pulse specifies the acoustic wave form. %%pulse = 1, 2, 3 and 4 correspond to sinusoidal wave, Gaussian wave, square wave and broadband wave, respectively. The implementation of the broadband wave is based on \cite Tam05

• %%npulse specifies the number of cycles of the acoustic wave generated. Only applies to %%pulse = 1 and 3 (sine and square waves), and must be an integer for non-planar waves.

• %%mag specifies the peak amplitude of the acoustic wave.

• %%frequency and %%wavelength specify the frequency and wavelength of the acoustic wave, respectively. These parameters are exclusive and exactly one of them must be specified for %%pulse = 1 or 3 (sine or square waves). They are related by the speed of sound in the medium: frequency = speed_of_sound / wavelength.

• %%gauss_sigma_time and %%gauss_sigma_dist specify the time and spatial widths of the Gaussian pulse in terms of sigma, respectively. In particular, %%gauss_sigma_time is the standard deviation in the Gaussian equation. These parameters are exclusive and exactly one of them must be specified for %%pulse = 2 (Gaussian wave). They are related by the speed of sound in the medium: gauss_sigma_dist = speed_of_sound * gauss_sigma_time.

• %%delay specifies the time delay of the acoustic wave. This parameter is optional for %%pulse = 1 or 3 (sine or square waves) and defaults to 0. It must be specified for %%pulse = 2 (Gaussian wave). It is important to note that setting the delay to 0 for a Gaussian pulse results in a half-Gaussian pulse, and delays that are too small may result in the pulse being cut off at the start of the simulation. 4*gauss_sigma_time is a typical value for the delay of a Gaussian pulse.

• %%dir specifies the direction of acoustic wave propagation for planar waves. The direction is defined by the angle in degrees from the x-axis in the x-y plane. It applies to both 2D and 3D simulation of planar waves (support is infinite in z-direction for 3D).

• %%length specifies the spatial length of the 2D or 3D planar wave. It is the length of the source plane perpendicular to the direction of wave propagation.

• %%height specifies the spatial height of the planar wave. Since %%dir is in the x-y plane, the height is perpendicular to the direction of wave propagation.

• %%foc_length specifies the focal length of the transducer for transducer waves. It is the distance from the transducer to the focal point.

• %%aperture specifies the aperture of the transducer. It is the diameter of the projection of the transducer arc onto the y-axis (2D) or spherical cap onto the y-z plane (3D). Set the aperture to double the focal length to simulate a transducer enclosing half of the circle/sphere. For the transducer array, it is the total aperture of the array.

• %%num_elements specifies the number of transducer elements in a transducer array.

• %%element_on specifies the element number of the transducer array that is on. The element number starts from 1, if all elements are on, set %%element_on to 0.

• %%element_spacing_angle specifies the spacing angle between adjacent transducers in radians. The total aperture (%%aperture) is set, so each transducer element is smaller if %%element_spacing_angle is larger.

• %%element_polygon_ratio specifies the ratio of the polygon side length to the aperture diameter of each transducer element in a circular 3D transducer array. The polygon side length is calculated by using the total aperture (%%aperture) as the circumcircle diameter and %%num_elements as the number of sides of the polygon. The ratio is used to specify the aperture size of each transducer element in the array as a ratio of the total aperture.

• %%rotate_angle specifies the rotation angle of the 3D circular transducer array along the x-axis (principal axis). It is optional and defaults to 0.

• %%bb_num_freq specifies the number discretized frequencies in the broadband acoustic wave. If %%bb_num_freq is 1, the acoustic wave will be a discrete tone (i.e. single frequency sine wave).

• %%bb_bandwidth specifies the bandwidth of the discretized frequencies.

• %%bb_lowest_freq specifies the lower frequency bound of the broadband acoustic wave. The upper frequency bound will be calculated as %%bb_lowest_freq + %%bb_num_freq * %%bb_bandwidth. The wave is no longer broadband below the lower bound and above the upper bound.

9. Sub-grid Bubble Models {#sec-bubble-models}

Parameter	Type	Description
bubbles_euler	Logical	Ensemble-averaged bubble modeling
bubbles_lagrange	Logical	Volume-averaged bubble modeling
bubble_model	Integer	[1] Gilmore; [2] Keller--Miksis; [3] Rayleigh-Plesset
Ca	Real	Cavitation number
Web	Real	Weber number
Re_inv	Real	Inverse Reynolds number
pref	Real	Reference pressure for bubble models
rhoref	Real	Reference density for bubble models
fluid_rho	Real	Reference fluid density
bub_pp%%R0ref*†‡	Real	Reference bubble radius
bub_pp%%p0ref*†‡	Real	Reference pressure
bub_pp%%rho0ref*†‡	Real	Reference density
bub_pp%%T0ref†‡	Real	Reference temperature
bub_pp%%ss*†‡	Real	Surface tension between host and gas (bubble)
bub_pp%%pv*†‡	Real	Vapor pressure of host
bub_pp%%vd† ‡	Real	Vapor diffusivity in gas (bubble)
bub_pp%%mu_l*†‡	Real	Viscosity of host in liquid state
bub_pp%%mu_v†	Real	Viscosity of host in vapor state
bub_pp%%mu_g†	Real	Viscosity of gas (bubble)
bub_pp%%gam_v†‡	Real	Specific heat ratio of host in vapor state
bub_pp%%gam_g*†‡	Real	Specific heat ratio of gas (bubble)
bub_pp%%M_v†‡	Real	Molecular weight of host
bub_pp%%M_g†‡	Real	Molecular weight of gas (bubble)
bub_pp%%k_v†‡	Real	Thermal conductivity of host in vapor state
bub_pp%%k_g†‡	Real	Thermal conductivity of gas (bubble)
bub_pp%%cp_v‡	Real	Specific heat capacity in constant pressure of host in vapor state
bub_pp%%cp_g‡	Real	Specific heat capacity in constant pressure of gas (bubble)
bub_pp%%R_v†‡	Real	Gas constant of host in vapor state
bub_pp%%R_g†‡	Real	Gas constant of gas (bubble)

• * These parameters are required for ensemble-averaged bubble model.
• † These parameters are required for ensemble-averaged bubble model with polytropic = F.
• ‡ These parameters are required for volume-averaged bubble model with polytropic = F.

This table lists the sub-grid bubble model parameters, which can be utilized in both the ensemble-averaged and volume-averaged bubble models.

• bubbles_euler activates the ensemble-averaged bubble model.

• bubbles_lagrange activates the volume-averaged bubble model.

• bubble_model specifies a model for spherical bubble dynamics by an integer of 1, 2, or 3. bubble_model = 1, 2, and 3 correspond to the Gilmore, Keller-Miksis, and Rayleigh-Plesset models.

• bub_pp specifies simulation parameters for the EE and/or EL bubble model.

Implementation of the parameters into the model follows \cite Ando10.

See @ref equations "Equations" Section 9 for the bubble dynamics equations.

9.1 Ensemble-Averaged Bubble Model

Parameter	Type	Description
bubbles_euler	Logical	Ensemble-averaged bubble modeling
polytropic	Logical	Polytropic gas compression
thermal	Integer	Thermal model: [1] Adiabatic; [2] Isothermal; [3] Transfer
polydisperse	Logical	Polydispersity in equilibrium bubble radius R0
nb	Integer	Number of bins: [1] Monodisperse; [$>1$] Polydisperse
poly_sigma	Real	Standard deviation for probability density function of polydisperse bubble populations
qbmm	Logical	Quadrature by method of moments
dist_type	Integer	Joint probability density function for bubble radius and velocity (only for qbmm = 'T')
sigR	Real	Standard deviation for the probability density function of bubble radius (only for qbmm = 'T') (EE)
sigV	Real	Standard deviation for the probability density function of bubble velocity (only for qbmm = 'T') (EE)
rhoRV	Real	Correlation coefficient for the joint probability density function of bubble radius and velocity (only for qbmm = 'T')

This table lists the ensemble-averaged bubble model parameters.

• polytropic activates polytropic gas compression in the bubble. When polytropic = 'F', the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on (\cite Preston07).

• thermal specifies a model for heat transfer across the bubble interface by an integer from 1 through 3. thermal = 1, 2, and 3 correspond to no heat transfer (adiabatic gas compression), isothermal heat transfer, and heat transfer with a constant heat transfer coefficient based on \cite Preston07, respectively.

• polydisperse activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius. Simpson's rule is used for integrating the log-normal PDF of equilibrium bubble radius for polydisperse populations.

• nb specifies the number of discrete bins that define the probability density function (PDF) of the equilibrium bubble radius.

• poly_sigma specifies the standard deviation of the log-normal PDF of equilibrium bubble radius for polydisperse populations.

• qbmm activates quadrature by method of moments, which assumes a PDF for bubble radius and velocity.

• dist_type specifies the initial joint PDF of initial bubble radius and bubble velocity required in qbmm. dist_type = 1 and 2 correspond to binormal and lognormal-normal distributions respectively.

• sigR specifies the standard deviation of the PDF of bubble radius required in the QBMM feature.

• sigV specifies the standard deviation of the PDF of bubble velocity required in the QBMM feature.

• rhoRV specifies the correlation coefficient of the joint PDF of bubble radius and bubble velocity required in the QBMM feature.

9.2 Volume-Averaged Bubble Model

Parameter	Type	Description
bubbles_lagrange	Logical	Lagrangian subgrid bubble model switch
nBubs_glb	Integer	Global number of bubbles
solver_approach	Integer	1: One-way coupling, 2: two-way coupling
cluster_type	Integer	Method to find p_inf
pressure_corrector	Logical	Cell pressure correction term
smooth_type	Integer	Smoothing function. 1: Gaussian, 2:Delta 3x3
heatTransfer_model	Logical	Activates the interface heat transfer model
massTransfer_model	Logical	Activates the interface mass transfer model
write_bubbles	Logical	Write files to track the bubble evolution each time step
write_bubbles_stats	Logical	Write the maximum and minimum radius of each bubble
epsilonb	Real	Standard deviation scaling for the gaussian function
charwidth	Real	Domain virtual depth (z direction, for 2D simulations)
valmaxvoid	Real	Maximum void fraction permitted

• nBubs_glb Total number of bubbles. Their initial conditions need to be specified in the ./input/lag_bubbles.dat file. See the example cases for additional information.

• solver_approach Specifies the Euler-Lagrange coupling method: [1] enables a one-way coupling approach, where the bubbles do not influence the Eulerian field. [2] activates the two-way coupling approach based on \cite Maeda18, where the effect of the bubbles is added in the Eulerian field as source terms.

• cluster_type Specifies method to find p_inf (pressure that drives the bubble dynamics): [1] activates the bilinear interpolation of the pressure field, while [2] enables the bubble dynamic closure based on \cite Maeda18, the full model is obtained when pressure_corrector is true.

• smooth_type Specifies the smoothening method of projecting the lagrangian bubbles in the Eulerian field: [1] activates the gaussian kernel function described in \cite Maeda18, while [2] activates the delta kernel function where the effect of the bubble is only seen in the specific bubble location cell.

• heatTransfer_model Activates the heat transfer model at the bubble's interface based on (\cite Preston07).

• massTransfer_model Activates the mass transfer model at the bubble's interface based on (\cite Preston07).

10. Velocity Field Setup {#sec-velocity-field-setup}

Parameter	Type	Description
perturb_flow	Logical	Perturb the initial velocity field by random noise
perturb_flow_fluid	Integer	Fluid density whose flow is to be perturbed
perturb_flow_mag	Real	Set the magnitude of flow perturbations
perturb_sph	Logical	Perturb the initial partial density by random noise
perturb_sph_fluid	Integer	Fluid component whose partial density is to be perturbed
mixlayer_vel_profile	Logical	Set the mean streamwise velocity to hyperbolic tangent profile
mixlayer_vel_coef	Real	Coefficient for the hyperbolic tangent profile of a mixing layer
mixlayer_perturb	Logical	Perturb the initial velocity field using a spectrum-based synthetic turbulence generation method
mixlayer_perturb_k0	Real	Base wavenumber for mixing layer perturbation
mixlayer_perturb_nk	Integer	Number of perturbation modes for mixing layer
simplex_perturb	Logical	Enable simplex noise perturbation of initial conditions

The table lists velocity field parameters. The parameters are optionally used to define initial velocity profiles and perturbations.

• perturb_flow activates the perturbation of initial velocity by random noise.

• perturb_flow_fluid specifies the fluid component whose flow will be perturbed.

• perturb_flow activates the perturbation of initial velocity by random noise.

• perturb_sph activates the perturbation of initial partial density by random noise.

• perturb_sph_fluid specifies the fluid component whose partial density is to be perturbed.

• mixlayer_vel_profile activates setting the mean streamwise velocity to a hyperbolic tangent profile. This option works only for n > 0.

• mixlayer_vel_coef is a parameter for the hyperbolic tangent profile of a mixing layer when mixlayer_vel_profile = 'T'. The mean streamwise velocity profile is given as:

\f[ u = \text{patch_icpp(1)%vel(1)} \cdot \tanh( y_{cc} \cdot \text{mixlayer_vel_coef}) \f]

• mixlayer_perturb activates the velocity perturbation for a temporal mixing layer with hyperbolic tangent mean streamwise velocity profile, using an inverter version of the spectrum-based synthetic turbulence generation method proposed by \cite Guo23. This option only works for p > 0 and mixlayer_vel_profile = 'T'.

11. Phase Change Model {#sec-phase-change}

Parameter	Type	Description
relax	Logical	Activates Phase Change model
relax_model	Integer	Phase change model: [5] pT-equilibrium; [6] pTg-equilibrium
palpha_eps	Real	tolerance of the Newton Solver to activate pT-equilibrium
ptgalpha_eps	Real	tolerance of the Newton Solver to activate pTg-equilibrium

• relax Activates the Phase Change model.

• relax_model Specifies the phase change model to be used: [5] enables pT-equilibrium, while [6] activates pTg-equilibrium (if criteria are met).

• palpha_eps Specifies the tolerance for the Newton Solvers used in the pT-equilibrium model.

• ptgalpha_eps Specifies the tolerance used for the Newton Solvers used in the pTg-equilibrium model.

12. Artificial Mach Number {#sec-artificial-mach-number}

Parameter	Type	Description
pi_fac	Real	Ratio of artificial and true pi_\infty values

• pi_fac specifies the ratio of artificial and true pi_\infty values (= artificial pi_\infty / true pi_\infty). This parameter enables the use of true pi_\infty in bubble dynamics models when the pi_\infty given in the case.py file is an artificial value.

13. Body Forces

Parameter	Type	Description
bf_x[y,z]	Logical	Enable body forces in the x[y,z] direction
k_x[y,y]	Real	Magnitude of oscillating acceleration
w_x[y,z]	Real	Frequency of oscillating acceleration
p_x[y,z]	Real	Phase shift of oscillating acceleration
g_x[y,z]	Real	Magnitude of background acceleration

k_x[y,z], w_x[y,z], p_x[y,z], and g_x[y,z] define an oscillating acceleration in the x[y,z] direction with the form

\f[ a_{x[y,z]} = g_{x[y,z]} + k_{x[y,z]}\sin\left(w_{x[y,z]}t + p_{x[y,z]}\right). \f]

By convention, positive accelerations in the x[y,z] direction are in the positive x[y,z] direction.

14. Magnetohydrodynamics (MHD)

Parameter	Type	Description
mhd	Logical	Enable ideal MHD simulation
relativity	Logical	Enable relativistic MHD simulation
hyper_cleaning	Logical	Enable hyperbolic (GLM) divergence cleaning for div B
hyper_cleaning_speed	Real	Cleaning wave speed c_h
hyper_cleaning_tau	Real	Cleaning damping timescale tau
Bx[y,z]	Real	Initial magnetic field in the x[y,z] direction
Bx0	Real	Constant magnetic field in the x direction (1D only)

• mhd is currently only available for single-component flows and 5-equation model. Its compatibility with most other features is work in progress.

• relativity only works for mhd enabled and activates relativistic MHD (RMHD) simulation.

• hyper_cleaning \cite Dedner02 only works with mhd in 2D/3D and reduces numerical div B errors by propagation and damping. Currently not compatible with HLLD (riemann_solver = 4).

• hyper_cleaning_speed sets the propagation speed of divergence-cleaning waves.

• hyper_cleaning_tau sets the decay timescale for divergence-cleaning.

• Bx0 is only used in 1D simulations to specify the constant magnetic field in the x direction. It must be specified in 1D simulations. Bx must not be used in 1D simulations.

• Bx, By, and Bz are used to specify the initial magnetic field in the x, y, and z directions, respectively. They must be specified in all 1D/2D/3D MHD simulations, with the exception of Bx in 1D simulations.

Note: In 1D/2D/3D simulations, all three velocity components are treated as state variables and must be specified in the case file.

Note: For relativistic flow, the conservative and primitive densities are different. rho_wrt outputs the primitive (rest mass) density.

15. Elasticity

Parameter	Type	Description
hypoelasticity	Logical	Enable hypoelasticity simulation
hyperelasticity	Logical	Enable hyperelasticity simulation
cont_damage	Logical	Enable continuum damage model
tau_star	Real	Threshold stress for continuum damage model
cont_damage_s	Real	Power s for continuum damage model
alpha_bar	Real	Damage factor (rate) for continuum damage model

• cont_damage activates continuum damage model for solid materials. Requires tau_star, cont_damage_s, and alpha_bar to be set (empirically determined) (\cite Cao19).

16. Cylindrical Coordinates

When cyl_coord = 'T' is set in 3D the following constraints must be met:

• bc_y%%beg = -14 enables the axis boundary condition

• bc_z%%beg = bc_z%%end = -1 enables periodic boundary conditions in the azimuthal direction

• z_domain%%beg = 0 sets the azimuthal starting point to 0

• z_domain%%end = 2*math.pi to set the azimuthal ending point to \f$2\pi\f$ (note, requires import math in the case file)

When cyl_coord = 'T' is set in 2D the following constraints must be met:

• bc_y%%beg = -2 to enable reflective boundary conditions

17. Chemistry

Parameter	Type	Description
chemistry	Logical	Enable chemistry simulation
chem_params%%diffusion	Logical	Enable multispecies diffusion
chem_params%%reactions	Logical	Enable chemical reactions
chem_params%%gamma_method	Integer	Methodology for calculating the heat capacity ratio
chem_params%%transport_model	Integer	Methodology for calculating the diffusion coefficients
cantera_file	String	Cantera-format mechanism file (e.g., .yaml)

• chem_params%%transport_model specifies the methodology for calculating diffusion coefficients and other transport properties, 1 for mixture-average, 2 for Unity-Lewis

• cantera_file specifies the chemical mechanism file. If the file is part of the standard Cantera library, only the filename is required. Otherwise, the file must be located in the same directory as your case.py file

18. Chemistry-Specific Boundary Conditions

Parameter	Type	Description
bc_[x,y,z]%%isothermal_in	Logical	Enable isothermal wall at the domain entrance (minimum coordinate).
bc_[x,y,z]%%isothermal_out	Logical	Enable isothermal wall at the domain exit (maximum coordinate).
bc_[x,y,z]%%Twall_in	Real	Temperature [K] of the entrance isothermal wall.
bc_[x,y,z]%%Twall_out	Real	Temperature [K] of the exit isothermal wall.

This boundary condition can be used for fixed-temperature (isothermal) walls at the domain extremities. It is exclusively available for reacting flows and requires chemistry to be enabled. It properly evaluates heat and species fluxes at the interface when chemistry = 'T', chem_params%%diffusion = 'T', and the corresponding domain boundary is set to a slip wall (bc_[x,y,z]%%[beg,end] = -15) or a no-slip wall (bc_[x,y,z]%%[beg,end] = -16).

19. Non-Newtonian (Herschel-Bulkley) Viscosity {#sec-non-newtonian}

Parameter	Type	Description
fluid_pp(i)%%non_newtonian	Logical	Enable Herschel-Bulkley non-Newtonian viscosity for fluid \f$i\f$.
fluid_pp(i)%%K	Real	Consistency index \f$K\f$.
fluid_pp(i)%%nn	Real	Flow index \f$n\f$ (\f$n<1\f$ shear-thinning, \f$n>1\f$ shear-thickening).
fluid_pp(i)%%tau0	Real	Yield stress \f$\tau_0\f$; set to 0 for pure power-law.
fluid_pp(i)%%hb_m	Real	Papanastasiou regularization parameter \f$m\f$; required when tau0 > 0.
fluid_pp(i)%%mu_min	Real	Lower viscosity clamp \f$\mu_{\min}\f$.
fluid_pp(i)%%mu_max	Real	Upper viscosity clamp \f$\mu_{\max}\f$ (required).
fluid_pp(i)%%mu_bulk	Real	Reserved; non-Newtonian bulk viscosity is not yet supported. The validator rejects it on a non-Newtonian fluid; on a Newtonian fluid it is accepted and ignored.

The effective dynamic viscosity is computed from the Papanastasiou-regularized Herschel-Bulkley model:

\f[ \mu_{\rm eff}(\dot\gamma) = \frac{\tau_0}{\dot\gamma}!\left(1 - e^{-m,\dot\gamma}\right) + K,\dot\gamma^{n-1}, \qquad \dot\gamma = \sqrt{2,D_{ij}D_{ij}}, \f]

where \f$D_{ij} = \frac{1}{2}(\partial_i u_j + \partial_j u_i)\f$ is the strain-rate tensor and \f$\dot\gamma\f$ is the scalar shear rate. The result is clamped to \f$[\mu_{\min},,\mu_{\max}]\f$.

Special cases:

• tau0 = 0: pure power-law fluid, \f$\mu_{\rm eff} = K,\dot\gamma^{n-1}\f$.
• tau0 = 0, nn = 1: Newtonian fluid with constant viscosity \f$\mu = K\f$.
• tau0 > 0, nn = 1: Bingham plastic.

Usage notes:

• Requires viscous = T. fluid_pp(i)%%Re(1) must be set (use 1.0/K to register the fluid as viscous; the HB model overrides \f$\mu_{\rm eff}\f$ cell-by-cell). fluid_pp(i)%%Re(2) (bulk viscosity) must not be set for a non-Newtonian fluid.
• mu_max is required; mu_min is inactive if omitted (no lower clamp applied).
• Positivity requirements: K, nn, and mu_max must be positive; mu_min and hb_m must be positive when set; tau0 must be non-negative.
• Requires model_eqns = 2 or 3 and is incompatible with igr.
• Supported only with riemann_solver = 1 (HLL) or riemann_solver = 2 (HLLC).
• The HB parameters (K, nn, tau0, hb_m, mu_min, mu_max, mu_bulk) may only be set on a fluid with non_newtonian = T; the validator rejects them otherwise.
• All HB parameters are non-dimensional (scaled by \f$\rho_{\rm ref} U_{\rm ref} L_{\rm ref}\f$), so \f$1/\mu_{\rm eff}\f$ equals the local effective Reynolds number.
• For cylindrical geometry (cyl_coord = T) the shear rate uses the grid-direction strain components; curvature corrections to \f$\dot\gamma\f$ are not yet included.

20. GPU Performance (NVIDIA UVM)

Parameter	Type	Description
nv_uvm_out_of_core	Logical	Enable NVIDIA Unified Virtual Memory out-of-core execution
nv_uvm_pref_gpu	Logical	Prefer GPU placement for NVIDIA UVM allocations
nv_uvm_igr_temps_on_gpu	Integer	Store IGR temporaries on GPU with UVM

• These parameters are for NVIDIA Grace-Hopper and similar architectures with hardware-managed unified memory. They allow MFC to run problems larger than GPU memory by paging data between host and device.

Enumerations

Boundary conditions {#boundary-conditions}

#	Type	Description
-1	Normal	Periodic
-2	Normal	Reflective
-3	Normal	Ghost cell extrapolation
-4	Normal	Riemann extrapolation
-5	Characteristic	Slip wall
-6	Characteristic	Non-reflecting subsonic buffer
-7	Characteristic	Non-reflecting subsonic inflow
-8	Characteristic	Non-reflecting subsonic outflow
-9	Characteristic	Force-free subsonic outflow
-10	Characteristic	Constant pressure subsonic outflow
-11	Characteristic	Supersonic inflow
-12	Characteristic	Supersonic outflow
-14	Normal	Axis *
-15	Normal	Slip wall
-16	Normal	No-slip wall

*: This boundary condition is only used for bc_y%%beg when using cylindrical coordinates (cyl_coord = 'T' and 3D). For axisymmetric problems, use bc_y%%beg = -2 with cyl_coord = 'T' in 2D.

The boundary condition supported by the MFC are listed in table Boundary Conditions. Their number (#) corresponds to the input value in input.py labeled bc_[x,y,z]%[beg,end] (see table Simulation Algorithm Parameters). The entries labeled "Characteristic." are characteristic boundary conditions based on \cite Thompson87 and \cite Thompson90.

Generalized Characteristic Boundary conditions

Parameter	Type	Description
bc_[x,y,z]%%grcbc_in	Logical	Enable grcbc for subsonic inflow
bc_[x,y,z]%%grcbc_out	Logical	Enable grcbc for subsonic outflow (pressure)
bc_[x,y,z]%%grcbc_vel_out	Logical	Enable grcbc for subsonic outflow (pressure + normal velocity)
bc_[x,y,z]%%vel_in	Real Array	Inflow velocities in x, y and z directions
bc_[x,y,z]%%vel_out	Real Array	Outflow velocities in x, y and z directions
bc_[x,y,z]%%pres_in	Real	Inflow pressure
bc_[x,y,z]%%pres_out	Real	Outflow pressure
bc_[x,y,z]%%alpha_rho_in	Real Array	Inflow density
bc_[x,y,z]%%alpha_in	Real Array	Inflow void fraction

This boundary condition can be used for subsonic inflow (bc_[x,y,z]%[beg,end] = -7) and subsonic outflow (bc_[x,y,z]%[beg,end] = -8) characteristic boundary conditions. These are based on \cite Pirozzoli13. This enables to provide inflow and outflow conditions outside the computational domain.

Patch types {#patch-types}

#	Name	Dim.	Smooth	Description
1	Line segment	1	N	Requires x_centroid and length_x.
2	Circle	2	Y	Requires [x,y]_centroid and radius.
3	Rectangle	2	N	Coordinate-aligned. Requires [x,y]_centroid and length_[x,y].
4	Sweep line	2	Y	Not coordinate aligned. Requires [x,y]_centroid and normal(i).
5	Ellipse	2	Y	Requires [x,y]_centroid and radii(i).
6	N/A	N/A	N/A	No longer exists. Empty.
7	N/A	N/A	N/A	No longer exists. Empty.
8	Sphere	3	Y	Requires [x,y,z]_centroid and radius
9	Cuboid	3	N	Coordinate-aligned. Requires [x,y,z]_centroid and length_[x,y,z].
10	Cylinder	3	Y	Requires [x,y,z]_centroid, radius, and length_[x,y,z].
11	Sweep plane	3	Y	Not coordinate-aligned. Requires x[y,z]_centroid and normal(i).
12	Ellipsoid	3	Y	Requires [x,y,z]_centroid and radii(i).
13	2D modal (Fourier)	2	Y	Requires x_centroid, y_centroid, radius. Optional: fourier_cos(n), fourier_sin(n) (n=1..10), modal_clip_r_to_min, modal_r_min, modal_use_exp_form.
14	3D spherical harmonic	3	Y	Requires x_centroid, y_centroid, z_centroid, radius. Optional: sph_har_coeff(l,m) (l=0..5, m=-l..l).
15	N/A	N/A	N/A	No longer exists. Empty.
16	1D bubble pulse	1	N	Requires x_centroid, length_x
17	Spiral	2	N	Requires [x,y]_centroid
18	2D Varcircle	2	Y	Requires [x,y]_centroid, radius, and thickness
19	3D Varcircle	3	Y	Requires [x,y,z]_centroid, length_z, radius, and thickness
20	2D Taylor-Green Vortex	2	N	Requires [x,y]_centroid, length_x, length_y, vel(1), and vel(2)
21	Model	2 & 3	Y	Imports a Model (STL/OBJ). Requires model_id.

The patch types supported by the MFC are listed in table Patch Types. This includes types exclusive to one-, two-, and three-dimensional problems. The patch type number (#) corresponds to the input value in input.py labeled patch_icpp(j)%%geometry where $j$ is the patch index. Each patch requires a different set of parameters, which are also listed in this table.

Geometry 13: 2D modal (Fourier):
Boundary is at polar angle \f$\theta = \mathrm{atan2}(y - y_{\mathrm{centroid}}, x - x_{\mathrm{centroid}})\f$.

• Additive form (default, modal_use_exp_form false):
  \f$R_{\mathrm{boundary}} = \mathrm{radius} + \sum_n \bigl[ \mathtt{fourier_cos}(n)\cos(n\theta) + \mathtt{fourier_sin}(n)\sin(n\theta) \bigr]\f$.
  Coefficients are absolute: same units as radius (length).
  If this formula gives \f$R_{\mathrm{boundary}} < 0\f$ at some \f$\theta\f$, it is clipped to zero.
  With modal_clip_r_to_min true, if \f$R_{\mathrm{boundary}} <\f$ modal_r_min at some \f$\theta\f$, it is clipped to modal_r_min.

• Exponential form (modal_use_exp_form true):
  \f$R_{\mathrm{boundary}} = \mathrm{radius} \times \exp\bigl( \sum_n [\ldots] \bigr)\f$.
  Coefficients are relative (dimensionless); the sum scales the radius.

Immersed Boundary Patch Types {#immersed-boundary-patch-types}

#	Name	Dim.
2	2D Circle	2
3	2D Rectangle	2
4	2D Airfoil	2
8	3D Sphere	3
10	3D Cylinder	3
11	3D Airfoil	3

Acoustic Supports {#acoustic-supports}

#	Name	Dim.	Requirements
1	Planar source	1D	%%loc(1), %%pulse, %%npulse, %%mag, and %%dir
2	Planar source	2D	#1 requirements, %%loc(2) and %%length
3	Planar source	3D	#2 requirements and %%height
5	Cylindrical Transducer	2D	%%loc(1), %%loc(2), %%pulse, %%npulse, %%mag, %%foc_length, %%aperture
6	Spherical Transducer	2D-Axisym	#5 requirements
7	Spherical Transducer	3D	#5 requirements and %%loc(3)
9	Arcuate Transducer Array	2D	#5 requirements, %%num_elements, %%element_on, %%element_spacing_angle
10	Annular Transducer Array	2D-Axisym	#9 requirements
11	Circular Transducer Array	3D	#7 requirements, %%element_polygon_ratio, and %%rotate_angle

The required parameters for each acoustic support type are listed in Acoustic Source. The acoustic support number (#) corresponds to the acoustic support type Acoustic(i)%%support, where $i$ is the acoustic source index. For each %%parameter, prepend the parameter with acoustic(i)%.

Additional requirements for all acoustic support types:

• acoustic_source = 'T' must be used to activate the acoustic source module.

• num_source must be set to the total number of acoustic sources.

• %%support must be set to the acoustic support number in the table.

• %%dipole is only supported for planar sources.

• %%npulse = 1 or 3 requires exactly one of %%frequency or %%wavelength to be set. It accepts %%delay as an optional parameter (default = 0).

• %%npulse = 2 requires exactly one of %%gauss_sigma_time or %%gauss_sigma_space to be set. It requires %%delay to be set.

Description of the acoustic support types:

• %%support = 1 specifies an infinite source plane that is normal to the $x$-axis and intersects with the axis at $x=$ %%loc(1) in 1D simulation. %%dir > 0 specifies a rightward propagating wave, and %%dir < 0 specifies a leftward propagating wave. %%dir = 0 is not allowed.

• %%support = 2 specifies a semi-infinite source plane in 2D simulation. The midplane location is [%%loc(1), %%loc(2)] and the normal vector is [\f$\mathrm{cos}\f$(%%dir), \f$\mathrm{sin}\f$(%%dir)]. The length of the source plane is %%length, and the plane is perpendicular to the direction of wave propagation (defined by %%dir).

• %%support = 3 specifies a semi-infinite source plane in 3D simulation. The midplane location is [%%loc(1), %%loc(2)] and the normal vector is [\f$\mathrm{cos}\f$(%%dir), \f$\mathrm{sin}\f$(%%dir)]. The length of the source plane is %%length, and the plane is perpendicular to the direction of wave propagation (defined by %%dir). Note that the plane is infinite in the $z$-direction, so %%loc(3) is not required.

• %%support = 5 specifies a circular transducer in 2D simulation. The transducer is centered at [%%loc(1), %%loc(2)] with a focal length of %%foc_length and an aperture of %%aperture. The center location is not the focal point; it is the tip of the circular arc (intersection of the arc and the x-axis). The aperture is the length of the projection of the circular arc onto the y-axis. If a semi-circle is desired, set the aperture to double the focal length. Note that this is physically a cylindrical transducer, and due to the complexity of Green's function for 2D wave, no closed-form solution is available for the 2D circular transducer, and an approximate is used (see \cite Maeda17 for details). For the mass source term correction factor, the theoretical approximation factor of -0.5 in ($r_{foc}^{-0.5}$) is replaced by an empirically determined factor of -0.85.

• %%support = 6 specifies a spherical transducer in 2D axisymmetric simulation. It is identical to %%support = 5 in terms of simulation parameters. Note that this is physically a spherical 3D transducer, so the equation is exact.

• %%support = 7 specifies a spherical transducer in 3D simulation. The transducer is centered at [%%loc(1), %%loc(2), %%loc(3)] with a focal length of %%foc_length and an aperture of %%aperture. The center location is not the focal point; it is the tip of the spherical cap (intersection of the cap and the x-axis). The aperture is the diameter of the projection of the spherical cap onto the y-z plane. If a semi-sphere is desired, set the aperture to double the focal length. Again, the equation is exact.

• %%support = 9 specifies an arcuate transducer array in 2D simulation. The total aperture of the array is %%aperture, which is similar to %%support = 5. The parameters %%num_elements and %%element_spacing_angle specify the number of transducer elements and the spacing angle. The spacing angle is the angle of the gap between adjacent transducer elements in the array. Because the total aperture is set, each transducer element is smaller if the spacing angle is larger. Physically it represents curved panels. Note that similar to %%support = 5, the mass source term correction factor is empirically determined to be -0.85.

• %%support = 10 specifies an annular transducer array in 2D axisymmetric simulation. It is identical to %%support = 9 in terms of simulation parameters. It physically represents the a annulus obtained by revolving the arc in %%support = 9 around the x-axis.

• %%support = 11 specifies a circular transducer array in 3D simulation. The total aperture of the array is %%aperture, which is similar to %%support = 7. The parameters %%num_elements, %%element_polygon_ratio, and %%rotate_angle specify the number of transducer elements, the ratio of the polygon side length to the transducer element radius, and the rotation angle of the array. The polygon side length is calculated by using the total aperture as the circumcicle diameter, and the number of sides of the polygon as %%num_elements. The ratio is used specify the aperture size of each transducer element in the array, as a ratio of the total aperture. The rotation angle is optional and defaults to 0. Physically it represents a circular ring of transducer elements.

Conservative Variables Ordering

5-eqn	6-eqn
num_fluids continuity variables	num_fluids continuity variables
num_dims momentum variables	num_dims momentum variables
1 energy variable	1 energy variable
num_fluids advection variables	num_fluids advection variables
N/A	num_fluids internal energy variables

The above variables are used for all simulations.

5-eqn	6-eqn
sub-grid bubble variables	N/A
hypoelastic variables	N/A

The above variables correspond to optional physics.

Primitive Variables Ordering

5-eqn	6-eqn
num_fluids densities	num_fluids densities
num_dims velocities	num_dims velocities
1 pressure	1 pressure
num_fluids volume fractions	num_fluids volume fractions
N/A	num_fluids partial pressures

The above variables are used for all simulations.

5-eqn	6-eqn
sub-grid bubble variables	N/A
hypoelastic variables	N/A

The above variables correspond to optional physics.

Page last updated: 2026-02-15