(using)=
# Using Ratel

The Ratel library includes support for processing inputs and handling command-line options.
Most examples and applications using Ratel will inherit these options.

## Command Line Options

Ratel is controlled via command-line options.
These command line options may be stored in a YML file specified by the runtime option `-options_file`.
The following command line options are mandatory:

:::{list-table} Mandatory Runtime Options
:header-rows: 1
:widths: 3 7

* - Option
  - Description

* - `-dm_plex_filename [filename]`
  - Path to mesh file in any format supported by PETSc.
    Alternatively, a built-in mesh, such as `-dm_plex_box` may be specified.

* - `-bc_clamp [int list]`
  - List of face sets on which to displace by `-bc_clamp_[facenumber]_translate [x,y,z,...]` and/or `bc_clamp_[facenumber]_rotate [rx,ry,rz,c_0,c_1]`.
    Note: The default for a clamped face is zero displacement.
    All displacement is with respect to the initial configuration.
    At least one fixed boundary is required to prevent a singular system.
:::

:::{note}
This solver can use any mesh format that PETSc's `DMPlex` can read (Exodus, Gmsh, Med, etc.).
Our tests have primarily been using Exodus meshes created using [CUBIT]; sample meshes used for the example runs suggested here can be found in [this repository].
Note that many mesh formats require PETSc to be configured appropriately; e.g., `--download-exodusii` for Exodus support.
:::

Consider the specific example of the mesh seen below:

```{image} https://github.com/jeremylt/ceedSampleMeshes/raw/master/cylinderDiagram.png
```

With the sidesets defined in the figure, we provide here an example of a minimal set of command line options:

```console
$ ./bin/ratel-quasistatic -dm_plex_filename [.exo file] -order 4 -E 1e6 -nu 0.3 -bc_clamp 998,999 -bc_clamp_998_translate 0,-0.5,1
```

In this example, we set the left boundary, face set $999$, to zero displacement and the right boundary, face set $998$, to displace $0$ in the $x$ direction, $-0.5$ in the $y$, and $1$ in the $z$.

As an alternative to specifying a mesh with {code}`-dm_plex_filename`, the user may use a DMPlex box mesh by specifying {code}`-dm_plex_box_faces [int list]`, {code}`-dm_plex_box_upper [real list]`, and {code}`-dm_plex_box_lower [real list]`.

As an alternative example exploiting {code}`-dm_plex_box_faces`, we consider a {code}`4 x 4 x 4` mesh where essential (Drichlet) boundary condition is placed on the top and bottom.
Side 1 is held in place while side 2 is rotated around $x$-axis:

```console
$ ./bin/ratel-quasistatic -model elasticity-neo-hookean-initial -E 1 -nu 0.3 -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 4,4,4 -bc_clamp 1,2 -bc_clamp_2_rotate 0,0,1,0,.05
```

:::{note}
If the coordinates for a particular side of a mesh are zero along the axis of rotation, it may appear that particular side is clamped zero.
:::

On each boundary node, the rotation magnitude is computed: {code}`theta = (c_0 + c_1 * cx) * loadIncrement` where {code}`cx = kx * x + ky * y + kz * z`, with {code}`kx`, {code}`ky`, {code}`kz` are normalized values.

The command line options just shown are the minimum requirements to run the application, but additional options may also be set as follows

:::{list-table} Additional Runtime Options
  :header-rows: 1
  :widths: 25 50 25

* - Option
  - Description
  - Default value

* - `-ceed [string]`
  - CEED resource specifier
  - `/cpu/self`

* - `-options_file [filename]`
  - Filepath to yml file with runtime options for materials with one field
  - `/path/to/my.yml`

* - `-method [fem | mpm]`
  - Numerical method used to run simulations, either Finite Element Method (`fem`) or Material Point Method (`mpm`). See `-mpm_*` for MPM specific options.
  - `fem`

* - `mpm_num_points [int]`
  - Total number of material points, divided evenly between cells in the domain. Only used if `-mpm_num_points_per_cell` is not provided.
  - `1728`

* - `-mpm_num_points_per_cell [int]`
  - Number of material points per cell, rounded up to the nearest cube. If not provided, `-mpm_num_points` is used.
  -

* - `-mpm_point_location_type [gauss | uniform | cell_random]`
  - Initial position of material points per-cell.
    If `gauss`, uses the Gauss quadrature points.
    If `uniform`, points are uniformly spaced within the cell. The distance between points across cell boundaries is also uniform.
    If `cell_random`, points are placed randomly with each cell. Note that each cell will still contain the same number of points.
  - `uniform`

* - `-order [int]`
  - Polynomial order of solution basis functions
  - `3`

* - `-orders [int list]`
  - Polynomial order of solution basis functions for all fields
  - `3`

* - `-q_extra [int]`
  - Increased quadrature space order; final order given by `order[0] + q_extra`
  - `0`

* - `-diagnostic_order [int]`
  - Order for diagnostic values mesh
  - Same value as multigrid fine level order specified via `-order`

* - `-diagnostic_geometry_order [int]`
  - Geometry order for diagnostic values mesh
  -

* - `-continue_file [filename]`
  - Filepath to binary file holding restart information and vector
  -

* - `-model`
  - Material model to use (`elasticity-linear`, `elasticity-neo-hookean-current`, `elasticity-mooney-rivlin-initial`, etc.)
  - `elasticity-linear`

* - `-forcing`
  - Forcing term option (`none`, `body`, or `mms`)
  -

* - `-forcing_body_vector [ax,ay,az][,a2x,a2y,a2z,...]`
  - User-specified acceleration vectors for applying body forces. Scaled by the mass computed with material density `-{material name}_rho` and nodal volume. If more than one vector is specified, transition times must be provided as `-forcing_body_times [t1,t2,...]`.
  -

* - `-forcing_body_times [t1,t2,...]`
  - Transition times between each acceleration vector.
    Note: If the first specified time is after `0`, an acceleration vector of `0,0,0` will be added with time `0` in order to ramp body forces.
    By default, linearly applies acceleration from `t = 0` to `t = 1`
  -

* - `-forcing_body_interpolation`
  - Interpolation type between specified acceleration vectors (`none` or `linear`).
  -

* - `-ceed_field_[fieldnumber]_mat_type [MatType]`
  - Submatrix `MatType` for a `MATCEED` submatrix for a single solution field
  - `MATCEED`

* - `-ceed_[fieldname]_mat_type [MatType]`
  - Submatrix `MatType` for a `MATCEED` submatrix for a single solution field
  - `MATCEED`

* - `-pc_type`
  - PC type for linear solver
  - `pmg` (most problems), `gamg` (linear problems with single active field), or `jacobi` (problems with multiple active fields)

* - `-pmg_coarsening_type [logarithmic | uniform | user]`
  - P-Multigrid coarsening to use (`logarithmic`, `uniform`, or `user`)
  - `logarithmic`

* - `-pmg_level_orders [int list]`
  - Array of orders for each multigrid level, in ascending order; fine grid order specified by `-orders`.
    This is only used with `-pmg_coarsening_type` `user`.
  -

* - `-pmg_coarse_order [int]`
  - Polynomial order of coarse grid basis functions for materials with one field.
    This is only used with `-pmg_coarsening_type` of `logarithmic` or `uniform`.
  - `1`

* - `-surface_force_faces [int or string list]`
  - List of named faces or face set label values on which to compute surface force and centroids.
  -

* - `-surface_force_face_[facename]_label_value [int]`
  - Face set label value for a given `facename`. Optional only if `facename` is an integer corresponding to a face set label value.
  -

* - `-surface_force_face_[facename]_bounding_box [xmin,ymin,zmin,xmax,ymax,zmax]`
  - Optional bounding box specifying a sub-region of a face set label value on which to compute surface forces and centroids. Useful for meshes which lack labeled faces.
  -

* - `-q_extra_surface_force [int]`
  - Increased quadrature space order for cell-to-face surface force computation; final order given by `order[0] + q_extra_surface_force`
  - `1`

* - `-expected_strain_energy [real]`
  - Expected strain energy, for testing
  -

* - `-expected_max_displacement_[direction] [displacement]`
  - Expected max displacement in a particular direction (`x`, `y`, or `z`), for testing
  -

* - `-expected_surface_force_face_[facenumber]_[direction] [force]`
  - Expected surface force on face `facenumber` in a particular direction (`x`, `y`, or `z`), for testing
  -

* - `-expected_surface_force_face_[facenumber] [fx,fy,fz]`
  - Expected surface forces on face `facenumber`, overrides any `-expected_surface_force_face_[facenumber]_[direction]`, for testing
  -

* - `-expected_centroid_face_[facenumber]_[direction] [value]`
  - Expected centroid of face `facenumber` in a particular direction (`x`, `y`, or `z`), for testing
  -

* - `-expected_centroid_face_[facenumber] [cx,cy,cz]`
  - Expected centroid of face `facenumber`, overrides any `-expected_centroid_face_[facenumber]`, for testing
  -

* - `-view_diagnostic_quantities viewer:filename.extension`
  - Output final solution for viewing, ASCII format to STDOUT is used if no viewer is passed
  -

* - `-view_swarm ascii:[base_filename].xmf`
  - Output final MPM fields at material points for viewing, only `XDMF` format is currently supported. Only supported for `-method mpm`.
  -

* - `-ts_monitor_strain_energy`
  - Output computed strain energy on each time step
  -

* - `-ts_monitor_diagnostic_quantities viewer:filename.extension`
  - Output final solution for viewing, ASCII format to STDOUT is used if no viewer is passed
  -

* - `-skip_diagnostic_quantities`
  - Skip diagnostic quantities computation all together after the simulation completes.
  -

* - `-ts_monitor_surface_force [ascii:filename.csv]`
  - Output centroid displacements and accurate surface forces computed using the volumetric residual operator (reaction force) for faces given by `-surface_force_faces` on each time step.
    Only CSV or plain text output is supported.
    If no filename is provided, results are printed to the terminal (`stdout`).
  -

* - `-ts_monitor_surface_force_per_face [ascii:base_filename.csv]`
  - Like `-ts_monitor_surface_force`, except that a separate file is created for each face provided to `-surface_force_faces`.
    The created files are named `[base filename]_[face name].csv`.
  -

* - `-ts_monitor_surface_force_cell_to_face`
  - Output centroid displacements and approximate surface forces computed using cell-to-face quadrature for faces given by `-surface_force_faces` on each time step.
    This option is extremely inaccurate depending on the boundary conditions and geometry.
    Generally, you should use `-ts_monitor_surface_force` instead.
  -

* - `-ts_monitor_swarm ascii:[base_filename].xmf`
  - Monitor all swarm fields (potentially very large) at each timestep, saving an XDMF file `[base_filename].xmf` which can be read by Paraview and
    one binary file `[base_filename]_k.pbin` per timestep `k` which are automatically referenced and indexed by the XDMF file.
    Only supported for `-method mpm`.
  -

* - `-ts_monitor_swarm_solution ascii:[base_filename].xmf`
  - Monitor swarm solution field at each timestep, saving an XDMF file `[base_filename].xmf` which can be read by Paraview and
    one binary file `[base_filename]_k.pbin` per timestep `k` which are automatically referenced and indexed by the XDMF file.
    Additional fields to monitor can be specified via `-ts_monitor_swarm_fields`.
    Only supported for `-method mpm`.
  -

* - `-ts_monitor_swarm_fields [list of field names]`
  - Specifies additional swarm fields to monitor when using `-ts_monitor_swarm_solution`.
    Reliably, all models define fields `J` and `rho` which can be used to compute the current mass density in postprocessing.
    Only supported for `-method mpm`.
  -

* - `-ts_monitor_checkpoint [filename]`
  - Output binary file with solution checkpoints for continuation.
    Note: Binary viewer and extension are always used.
  -

* - `-ts_monitor_checkpoint_toggle`
  - Toggle between writing two checkpoint files, `filename-0.bin` and `filename-1.bin`.
  -

* - `-checkpoint_final_solution [filename]`
  - Output binary file with final solution checkpoint.
    Note: Automatic with `ts_monitor_checkpoint`.
  -

* - `-ts_monitor_*_interval [int]`
  - Number of time steps between monitor output, for any of `-ts_monitor_*` options above
  - `1`

* - `-ts_view`
  - View PETSc `TS` time stepper configuration
  -

* - `-snes_view`
  - View PETSc `SNES` nonlinear solver configuration
  -

* - `-log_view`
  - View PETSc performance log
  -

* - `-ts_monitor_reaction_force`
  - Deprecated: Use `-ts_monitor_surface_force`
  -

* - `-help`
  - View comprehensive information about run-time options
  -
:::

To verify the convergence of the linear elasticity formulation on a given mesh with the method of manufactured solutions, run:

```console
$ ./bin/ratel-quasistatic -dm_plex_filename [mesh] -order [order] -model elasticity-linear -nu [nu] -E [E] -forcing mms
```

This option attempts to recover a known solution from an analytically computed forcing term.

## Material Properties

Each material model has properties that need to be specified.
All properties are mandatory.

:::{list-table} Material Property and Model Options
:header-rows: 1

* - Option
  - Description
  - Model

* - `-E [real]`
  - [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus), $E > 0$
  - Neo-Hookean

* - `-nu [real]`
  - [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu < 0.5$
  - Neo-Hookean, Mooney-Rivlin or Ogden

* - `-nu_smoother [real]`
  - Poisson's ratio for multigrid smoothers, $\nu < 0.5$
  - Neo-Hookean, Mooney-Rivlin or Ogden

* - `-mu_1 [real]`
  - [Mooney-Rivlin material constant](https://en.wikipedia.org/wiki/Mooney–Rivlin_solid), $\mu_1 > 0$,
  - Mooney-Rivlin

* - `-mu_2 [real]`
  - [Mooney-Rivlin material constant](https://en.wikipedia.org/wiki/Mooney–Rivlin_solid), $\mu_2 > 0$
  - Mooney-Rivlin

* - `-m [m_1, m_2, m_3]`
  - [Ogden material constant](https://en.wikipedia.org/wiki/Ogden_hyperelastic_model)
  - Ogden

* - `-alpha [alpha_1, alpha_2, alpha_3]`
  - [Ogden material constant](https://en.wikipedia.org/wiki/Ogden_hyperelastic_model), $2 \secondlame = \sum_{j=1}^N m_j \alpha_j,$ with $m_j \alpha_j > 0$
  - Ogden

* - `deviatoric_viscosity [real]`
  - Deviatoric viscosity, $\eta_d > 0$
  - Viscoelasticity

* - `-yield_stress [real]`
  - [Initial yield stress threshold](https://en.wikipedia.org/wiki/Plasticity_(physics)), $\sigma_0 > 0$
  - Plasticity

* - `-linear_hardening [real]`
  - [Isotropic linear hardening parameter](https://en.wikipedia.org/wiki/Plasticity_(physics)), $H >= 0$
  - Plasticity

* - `-saturation_stress [real]`
  - [Saturation stress for Voce hardening component](https://en.wikipedia.org/wiki/Plasticity_(physics)), $\sigma_\infty >= \sigma_0$
  - Plasticity

* - `-hardening_decay [real]`
  - [Hardening decay parameter for Voce hardening component](https://en.wikipedia.org/wiki/Plasticity_(physics)), $\beta >= 0$
  - Plasticity

* - `-use_AT1 [bool]`
  - True (default): AT1 model, False: AT2
  - Phase Field Modeling of Fracture

* - `-use_offdiagonal [bool]`
  - True (default): full Jacobian, False: Jacobian diagonal terms only (for quasi-Newton methods),
  - Phase Field Modeling of Fracture

* - `-fracture_toughness [real]`
  - [Fracture toughness expressed as critical energy release rate](https://en.wikipedia.org/wiki/Energy_release_rate_(fracture_mechanics)), $G_c > 0$
  - Phase Field Modeling of Fracture

* - `-characteristic_length [real]`
  - Length scale parameter of phase-field model, $l_0 > 0$
  - Phase Field Modeling of Fracture

* - `-residual_stiffness [real]`
  - Residual stiffness in full fracture state ($\phi = 1$), $\eta > 0$
  - Phase Field Modeling of Fracture

* - `-damage_viscosity [real]`
  - Viscosity for viscous regularization of damage problem, $\zeta >= 0$
  - Phase Field Modeling of Fracture

* - `damage_scaling [real]`
  - Scaling parameter, $c_0 > 0$
  - Phase Field Modeling of Fracture

* - `-kappa [real]`
  - Fiber dispersion parameter, $\kappa\in[0,1/3]$. See {ref}`hgo-fibers` | Hyperbolic porosity model exponent, $\kappa >= 0$
  - Holzapfel-Gasser-Ogden | Poromechanics, finite strain theory

* - `-k1 [positive real]`
  - Stress-like anisotropic parameter, $k_1 > 0$
  - Holzapfel-Gasser-Ogden

* - `-k2 [positive real]`
  - Dimensionless exponential coeficient anisotropic parameter, $k_2 > 0$
  - Holzapfel-Gasser-Ogden

* - `-directions [x1,y1,z1][,x2,y2,z2][,x3,y3,z3] `
  - Reference directions for transversely anisotropic fiber bundles
  Up to 3 directions may be specified
  Directions do not need to be normalized or orthogonal, but should not be collinear
  - Holzapfel-Gasser-Ogden

* - `-lambda_d [real]`
  - First Lame parameter for drained mixture (solid skeleton), $\lambda_d >= 0$
  - Poromechanics

* - `-mu_d [real]`
  - Shear modulus for drained mixture (solid skeleton), $\mu_d >= 0$
  - Poromechanics

* - `-bulk_d [real]`
  - Bulk modulus for drained mixture (solid skeleton), $\bulk_d >= 0$
  - Poromechanics

* - `-bulk_s [real]`
  - Bulk modulus for solid, $\bulk_\rs >= 0$
  - Poromechanics

* - `-bulk_f [real]`
  - Bulk modulus for pore fluid, $\bulk_\rf >= 0$
  - Poromechanics

* - `-rho [real]`
  - Real solid mass density, $\rho^\sR_0 > 0$
  - Poromechanics

* - `-rho_fR0 [real]`
  - Real pore fluid mass density, $\rho^\fR_0 > 0$
  - Poromechanics

* - `-phi [real]`
  - Porosity, $0 < \phi^\rf < 1$
  - Poromechanics

* - `-mu_f [real]`
  - Pore fluid viscosity, $\mu_\rf >= 0$
  - Poromechanics

* - `-varkappa_0 [real]`
  - Intrinsic permeability, $\varkappa_{0} >= 0$
  - Poromechanics

* - `-pf_0 [real]`
  - Initial pore fluid pressure (required $> 0$ when $\Gamma > 0$)
  - Poromechanics

* - `-Gamma [real]`
  - Ideal gas ratio of specific heats, $\Gamma > 0$
  - Poromechanics, finite strain theory

* - `alpha_stab [real]`
  - Pore fluid pressure stabilization, $\alpha^\text{stab} > 0$
  - Poromechanics, finite strain theory

* - `-alpha [real]`
  - Biot's coefficient, $0 <= \alpha <= 1$
  - Poromechanics, linear theory

* - `-M [real]`
  - Biot's modulus, $M >= 0$
  - Poromechanics, linear theory
:::

(multiple-materials)=
## Multiple Materials

Ratel supports the use of solving with different material models defined for different segments of the mesh.
This feature requires additional runtime flags as well as some modifications to existing flags.
Different materials should be specified over labeled volumes of the mesh; an example of the header of a Gmsh mesh (provided in {code}`examples/meshes/materials_2.msh`) with two materials ("rod" and "binder") is shown below:
```
$ head examples/meshes/materials_2.msh
$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
4
2 1 "start"
2 2 "end"
3 3 "rod"
3 4 "binder"
$EndPhysicalNames
```
In this example, the ID value of the "rod" and "binder" volumes are 3 and 4, respectively.
In order to tell Ratel to treat these volumes as different materials, we use {code}`-material rod,binder` to provide label names for our specified materials (Note: these names do not have to match the names in the Gmsh mesh).
These label names will be used as prefixes (as {code}`-{material name}_`) to specify other aspects for each material at runtime.
We also specify, for each material, which domain label values to use with {code}`-rod_label_value 3 -binder_label_value 4`.
To define material parameters such as $E$ and $\nu$, we now use {code}`-binder_E 2.0 -binder_nu 0.4`.

An example set of command line options for the setting rods and binder materials is given below:

```console
$ ./bin/ratel-quasistatic -material rod,binder -rod_model elasticity-mooney-rivlin-initial -rod_mu_1 0.5 -rod_mu_2 0.5 -rod_nu 0.4 -binder_label_value 3 -binder_model elasticity-neo-hookean-initial -binder_E 2.0 -binder_nu 0.4 -binder_label_value 4
```

A complete list of command line options for specifying multiple materials is given below in the next table:

:::{list-table} Multiple Material Runtime Options
:header-rows: 1

* - Option
  - Description
  - Default value

* - `-material [string list]`
  - List of names to use as labels for each material.
  -

* - `-{material name}_model [model]`
  - Material to use (`elasticity-linear`, `elasticity-neo-hookean-current`, `elasticity-mooney-rivlin-initial`, etc.) for a given material.
  - `elasticity-linear`

* - `-{material name}_label_name [str]`
  - Domain label specifying the type of volume to use for specifying materials.
    Optional.
  - `"Cell Sets"`

* - `-{material name}_label_value [int]`
  - Domain value specifying the volume to use for a given material.
  -

* - `-{material name}_rho [real]`
  - Density of materal, by default in kg/m^3
  - `1`

* - `-{material name}_E [real]`
  - [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus), $E > 0$
  -

* - `-{material name}_nu [real]`
  - [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu < 0.5$
  -

* - `-{material name}_nu_smoother [real]`
  - Poisson's ratio for multigrid smoothers, $\nu < 0.5$
  -

* - `-{material name}_forcing [none | body | mms]`
  - Forcing term option (`none`, `body`, or `mms`).
    Overrides global forcing option provided with `-forcing`.
  - `none`

* - `-{material name}_forcing_body_vector [ax,ay,az][,a2x,a2y,a2z,...]`
  - User-specified acceleration vectors for applying body forces.
    Scaled by the mass computed with material density `-{material name}_rho` and nodal volume.
    If more than one vector is specified, transition times must be provided as `-forcing_body_times [t1,t2,...]`.
    Overrides global forcing option provided with `-forcing_body_vector`.
  - `0,0,0`

* - `-{material name}_forcing_body_times [t1,t2,...]`
  - Transition times between each acceleration vector.
    Note: If the first specified time is after `0`, an acceleration vector of `0,0,0` will be added with time `0` in order to ramp body forces.
    By default, linearly applies acceleration from `t = 0` to `t = 1`.
    Overrides global forcing option provided with `-forcing_body_times`.
  - `1`

* - `-{material name}_forcing_body_interpolation [none | linear]`
  - Interpolation type between specified acceleration vectors (`none` or `linear`).
    Overrides global forcing option provided with `-forcing_body_interpolation`.
  - `linear`
:::

An example of specifying a two material quasistatic solve with YAML is provided in {code}`examples/ex02-quasistatic-elasticity-multi-material.yml`.

## Algebraic Solvers

The examples are configured to use the following Newton-Krylov-Multigrid method by default.

- Newton-type methods for the nonlinear solve, with the hyperelasticity models globalized using load increments.
- Preconditioned conjugate gradients to solve the symmetric positive definite linear systems arising at each Newton step.
- Preconditioning via $p$-version multigrid coarsening to linear elements, with algebraic multigrid (PETSc's `GAMG`) for the coarse solve.
  The default smoother uses degree 3 Chebyshev with Jacobi preconditioning.
  (Lower degree is often faster, albeit less robust; try {code}`-mg_levels_ksp_max_it 2`, for example.)
  Application of the linear operators for all levels with order $p > 1$ is performed matrix-free using analytic Newton linearization, while the lowest order $p = 1$ operators are assembled explicitly (using coloring at present).

Many related solvers can be implemented by composing PETSc command-line options.
For example, to use Hypre's BoomerAMG for the coarse solve (using the assembled linear elements), one would use `-mg_coarse_pc_type hypre`.
Run with `-help` to see (many!) available command-line options related to algebraic solvers.

## Nondimensionalization

Quantities such as the Young's modulus vary over many orders of magnitude, and thus can lead to poorly scaled equations.
One can nondimensionalize the model by choosing an alternate system of units, such that displacements and residuals are of reasonable scales.

:::{list-table} (Non)dimensionalization options
:header-rows: 1

* - Option
  - Description
  - Default value

* - `-units_meter`
  - 1 meter in scaled length units
  - `1`

* - `-units_second`
  - 1 second in scaled time units
  - `1`

* - `-units_kilogram`
  - 1 kilogram in scaled mass units
  - `1`
:::

For example, consider a problem involving metals subject to gravity.

:::{list-table} Characteristic units for metals
:header-rows: 1

* - Quantity
  - Typical value in SI units

* - Displacement, $\bm u$
  - $1 \,\mathrm{cm} = 10^{-2} \,\mathrm m$

* - Young's modulus, $E$
  - $10^{11} \,\mathrm{Pa} = 10^{11} \,\mathrm{kg}\, \mathrm{m}^{-1}\, \mathrm s^{-2}$

* - Body force (gravity) on volume, $\int \rho \bm g$
  - $5 \cdot 10^4 \,\mathrm{kg}\, \mathrm m^{-2} \, \mathrm s^{-2} \cdot (\text{volume} \, \mathrm m^3)$
:::

One can choose units of displacement independently (e.g., {code}`-units_meter 100` to measure displacement in centimeters), but $E$ and $\int \rho \bm g$ have the same dependence on mass and time, so cannot both be made of order 1.
This reflects the fact that both quantities are not equally significant for a given displacement size; the relative significance of gravity increases as the domain size grows.

## Diagnostic Quantities

Diagnostic quantities for viewing are provided when the command line option for visualization output, {code}`-view_diagnostic_quantities viewer:filename.extension` is used.
The diagnostic quantities include displacement in the $x$ direction, displacement in the $y$ direction, displacement in the $z$ direction, pressure, $\trace \bm{E}$, $\trace \bm{E}^2$, $ J $, and strain energy density $\psi$.
The table below summarizes the formulations of each of these quantities for each material type.

:::{list-table} Diagnostic quantities
   :header-rows: 1

   * - Quantity
     - Linear elasticity
     - Neo-Hookean hyperelasticity

   * - Cauchy stress tensor
     - $\bm{\sigma}$
     - $\bm{\sigma} = \bm{F} \bm{S} \bm{F}^T / J $

   * - Hydrostatic pressure
     - $-\trace \bm{\sigma} / 3$
     - $-\trace \bm{\sigma} / 3$

   * - Volumetric strain
     - $\trace \bm{\varepsilon}$
     - $\trace \bm{E}$

   * - $\trace \bm{E}^2$
     - $\trace \bm{\varepsilon}^2$
     - $\trace \bm{E}^2$

   * - $ J $
     - $1 + \trace \bm{\varepsilon}$
     - $ J $

   * - Strain energy density
     - $\frac{\lambda}{2} (\trace \bm{\varepsilon})^2 + \mu \bm{\varepsilon} : \bm{\varepsilon}$
     - $\frac{\lambda}{4}\left(J^2 - 1 - 2 \log J \right) + \mu \trace \bm{E} - \mu \log J$

   * - von Mises stress
     - $\sqrt{\frac 3 2 \bm{\hat \sigma} \tcolon \bm{\hat \sigma}}$
     - $\sqrt{\frac 3 2 \bm{\hat \sigma} \tcolon \bm{\hat \sigma}}$

   * - Mass Density
     - $\rho_0 / J$
     - $\rho_0 / J$
:::

The von Mises stress uses the deviatoric part $\bm{\hat\sigma} = \bm{\sigma} - \frac 1 3 \trace \bm{\sigma}$ of the Cauchy stress $\bm{\sigma}$.
The integrated strain energy $\Psi = \int_{\Omega_0} \psi$ is also computed and printed upon completion of a solve.

Alternatively, one can skip all diagnostic evaluations with `-skip_diagnostic_quantities` if no diagnostic quantities are needed after the simulation finishes running.

[cubit]: https://cubit.sandia.gov/
[this repository]: https://github.com/jeremylt/ceedSampleMeshes

## Boundary Conditions

At least one Dirichlet/essential boundary condition is required in order to constrain the rigid body modes of the displacement field.
If you know that the rigid body modes are constrained by other means, for example by using frictional contact, use the flag `-bc_allow_no_clamp` to skip checking for essential boundaries

(using-dirichlet)=

### Dirichlet/Essential Boundaries
Ratel offers two main forms of Dirichlet or essential boundaries: `clamp` and `slip`.
To add `slip` or `clamp` boundary conditions on the first field of a model (displacement), use the options:
```
-bc_slip [list of facenumbers]
-bc_clamp [list of facenumbers]
```
To constrain a field other than displacement, only slip boundaries are supported and the field index or field name must also be specified via:
```
-bc_slip_[fieldname] [list of facenumbers]
-bc_slip_field_[fieldindex] [list of facenumbers]
```

#### Clamp

Clamp boundary conditions are a special displacement-only boundary condition which support prescribed translation and rotation of a given face.
If the first field has more than three components, for example a damage component, only the displacement degrees-of-freedom are constrained.
By default, clamp boundaries are homogenous Dirichlet boundaries, that is, the prescribed displacement is `0,0,0`.
The [Clamp Options Table](table_clamp_options) explains the command-line options available for clamp boundary conditions.
For brevity, the string `[prefix]` is used in place of the full prefix `bc_clamp_[facenumber]` in the options listings.

(table_clamp_options)=

:::{list-table} Clamp Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-bc_clamp [int list]`
  - List of face sets on which to apply clamp boundaries with prefix `bc_clamp_[facenumber]`.
    Note: The default for a clamped face is zero displacement.
  -

* - `-[prefix]_translate [x,y,z][,x2,y2,z2,...]`
  - One or more vectors specifying rigid translation of the face.
    Note: If more than one vector is specified, transition times must be provided as `-[prefix]_times [t1,t2,...]`.
    Note: If the constrained field has more than 3 components, such as a damage component, only the displacement components are constrained.
  - `0,0,0`

* - `-[prefix]_rotate [rx,ry,rz,c_0,c_1][,r2x,r2y,r2z,c2_0,c2_1,...]`
  - One or more rotation axes and rotation polynomial coefficients specifying rigid rotation of the face.
    Note: If more than one rotation is specified, transition times must be provided as `-[prefix]_times [t1,t2,...]`.
    Note: If the constrained field has more than 3 components, such as a damage component, only the displacement components are constrained.
  - `0,0,0,0,0`

* - `-[prefix]_times [t1,t2,...]`
  - Transition times between each rigid rotation/translation.
    Note: If the first specified time is after `0`, a rotation and translation of `0,0,0,0,0` and `0,0,0`, respectively, will be added with time `0` in order to start from an undeformed configuration.
  - `1`

* - `-[prefix]_interpolation [none | linear]`
  - Interpolation type between specified rotations/translations (`none` or `linear`).
  - `linear`
:::

#### Slip

Slip boundary conditions are simpler and more flexible than their clamp counterparts.
Each slip boundary requires that the constrained component indices (local to the constrained field) are listed explicitly.
As such, slip boundaries can be used to implement common engineering boundary conditions, such as rollers and symmetry boundaries.
As mentioned above, slip boundary option `[option]` for field `fieldname` with field index `fieldnumber` can be specified by:
```
-bc_slip_[fieldname]_face_[facenumber]_[option]
-bc_slip_field_[fieldnumber]_face_[facenumber]_[option]
```
For single-field models (or to only constrain the first field), the following simplified option can be used:
```
-bc_slip_[facenumber]_[option]
```
The [Slip Options Table](table_slip_options) explains the command-line options available for slip boundary conditions.
For brevity, the string `[prefix]` is used in place of one of the three prefixes above (`bc_slip_[facenumber]`, `bc_slip_field_[fieldnumber]_face_[facenumber]`, or `bc_slip_[fieldname]_face_[facenumber]`) in the options listings.

(table_slip_options)=

:::{list-table} Slip Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-bc_slip [int list]`
  - List of face sets on which to set slip boundary conditions with the prefix `-bc_slip_[facenumber]` on the field with index `0` (usually displacement or displacement+damage).
  -

* - `-bc_slip_field_[fieldnumber] [int list]`
  - List of face sets on which to set slip boundary conditions with the prefix `bc_slip_field_[fieldnumber]_face_[facenumber]` for solution field with index `fieldnumber`.
  -

* - `-bc_slip_[fieldname] [int list]`
  - List of face sets on which to set slip boundary conditions with the prefix `bc_slip_[fieldname]_face_[facenumber]` for solution field with name `fieldname`.
  -

* - `-[prefix]_components [int list]`
  - List of indices for the $k$ components to constrain on the face. The order in which component indices are provided is the same as the order of vector components for `-bc_slip_[facenumber]_translate`.
  -

* - `-[prefix]_translate [c1,...,ck][c1_2,...,ck_2,...,c1_n,...,ck_n]`
  - One or more vectors specifying prescribed translation of the $k$ constrained components of the face. That is, each translation vector should have length $k$.
    Note: If more than one vector is specified, transition times must be provided as `bc_slip_[facenumber]_times [t1,t2,...]`.
  - Zero vector of length $k$

* - `-[prefix]_times [t1,t2,...]`
  - Transition times between each prescribed translation.
    Note: If the first specified time is after `0`, a translation of `0,0,0` will be implicitly added with time `0` in order to start from an undeformed configuration.
  - `1`

* - `-[prefix]_interpolation [none | linear]`
  - Interpolation type between prescribed translations (`none` or `linear`).
  - `linear`
:::

#### Method of Manufactured Solutions (MMS)

Select material models support method of manufactured solutions problems for code verification.
Such MMS problems require external forcing terms, which can be applied via the `-bc_mms` interface.
Note: these are unlikely to be useful beyond specific test cases in the Ratel test suite.
The [MMS Options Table](table_mms_options) explains the command-line options available for MMS boundary conditions.

(table_mms_options)=

:::{list-table} MMS Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-bc_mms [int list]`
  - List of face sets on which to set Dirichlet boundary conditions to match a MMS solution
  -

* - `-bc_mms_field_[fieldnumber] [int list]`
  - List of face sets on which to set Dirichlet boundary conditions to match a MMS solution for a single solution field
  -

* - `-bc_mms_[fieldname] [int list]`
  - List of face sets on which to set Dirichlet boundary conditions to match a MMS solution for a single solution field
  -
:::

(using-traction)=

### Traction

Applies face traction vectors in the initial configuration can be specified via the `-bc_traction` interface.
These boundaries are the traditional Neumann boundary conditions with respect to the initial configuration.
The [Traction Options Table](table_traction_options) explains the command-line options available for contact boundary conditions.
For brevity, the string `[prefix]` is used in place of the full prefix `bc_traction_[facenumber]` in the option listings.


(table_traction_options)=

:::{list-table} Traction Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-bc_traction [int list]`
  - List of face sets on which to set traction boundary conditions with prefix `bc_traction_[facenumber]`.
  -

* - `-[prefix] [tx,ty,tz][,t2x,t2y,t2z,...]`
  - Traction vector(s) for face with components given with respect to the global coordinate system.
    If more than one vector is specified, transition times must be provided as `-[prefix]_times [t1,t2,...]`.
    Pressure-like units.
  -

* - `-[prefix]_times [t1,t2,...]`
  - Transition times between each traction vector.
    Note: If the first specified time is after `0`, a traction vector of `0,0,0` will be added with time `0` in order to start from an undeformed configuration.
    By default, linearly applies traction from `t = 0` to `t = 1`
  - `1`

* - `-[prefix]_interpolation [none | linear]`
  - Interpolation type between specified traction vectors (`none` or `linear`).
  - `linear`
:::

(using-pressure)=

### Pressure

Applied pressure boundary conditions can be specified via the `-bc_pressure` interface.
For details on the formulation of pressure boundaries, see {ref}`formulation_pressure_bc`.
The [Pressure Options Table](table_pressure_options) explains the command-line options available for contact boundary conditions.
For brevity, the string `[prefix]` is used in place of the full prefix `bc_pressure_[facenumber]` in the option listings.


(table_pressure_options)=

:::{list-table} Pressure Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-bc_pressure [int list]`
  - List of face sets on which to set pressure boundary conditions with prefix `bc_pressure_[facenumber]`.
  -

* - `-[prefix] [real t1p][,t2p,t3p,...]`
  - Pressure value(s) for face with components given with respect to the global coordinate system.
    If more than one pressure is specified, transition times must be provided as `-bc_pressure_[facenumber]_times [t1,t2,...]`.
  -

* - `-[prefix]_times [t1,t2,...]`
  - Transition times between each pressure value.
    Note: If the first specified time is after `0`, a pressure vector of `0` will be added with time `0` in order to start from an undeformed configuration.
    By default, linearly applies pressure from `t = 0` to `t = 1`
  - `1`

* - `-[prefix]_interpolation [none | linear]`
  - Interpolation type between specified pressure values (`none` or `linear`).
  - `linear`
:::

(using-contact)=

### Contact
Ratel supports frictionless and frictional contact with a rigid shapes via the `-bc_contact` command line option.
See {ref}`contact-boundary-conditions` for details on the theory of Ratel's contact enforcement methods and
see {ref}`friction-models` for details on the different friction models available in Ratel.

The [Contact Options Table](table_contact_options) explains the command-line options available for contact boundary conditions.
For brevity, the string `[prefix]` is used in place of the full prefix `bc_contact_[face_name]` in the option listings.

(table_contact_options)=

:::{list-table} Contact Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-bc_contact [string list]`
  - List of labels or names for each platen (half-plane) contact boundary condition to apply. Examples: "top,bottom" or "1,2".
    If integer values are provided, then they are assumed to correspond to surface attribute tags in the mesh.
    Otherwise, `-[prefix]_label_value [int]` must be set for each face with a non-integer name.
  -

* - `-[prefix]_shape [platen | cylinder]`
  - Specify the rigid shape of the contact surface, see {ref}`using-contact-shapes` for additional options for each type of shape.
    For historical reasons, `platen` is the default shape.
  - `platen`

* - `-[prefix]_type [nitsche | penalty]`
  - Solver method to use for platen boundary, either `nitsche` for Nitsche's method or `penalty` for penalty method.
  - `nitsche`

* - `-[prefix]_label_value [int]`
  - Value of `DMPlex` label to apply the platen to (such as the face label). If `[face_name]` is a string then the label value must be provided.
  - `[face_name]` if `[face_name]` is an integer.

* - `-[prefix]_translate [x,y,z][,x2,y2,z2,...]`
  - One or more vectors specifying rigid translation of the shape.
    Note: If more than one vector is specified, transition times must be provided as `[prefix]_times [t1,t2,...]`.
  - `0,0,0`

* - `-[prefix]_times [t1,t2,...]`
  - Transition times between each translation value, primarily used to control the shape's velocity.
    Note: If the first specified time is after `0`, a distance of `0` will be added implicitly with time `0` in order to start from an undeformed configuration.
    Note: If a single translation vector is provided, the time value defaults to `1`; i.e., the shape will move the prescribed translation in `1` time unit.
  - `1`

* - `-[prefix]_interpolation [none | linear]`
  - Interpolation type between specified distance values (`none` or `linear`).
    Note: `none` results in piece-wise constant displacement while `linear` results in piece-wise linear displacement and piece-wise constant velocity.
  - `linear`

* - `-[prefix]_gamma [real]`
  - Nitsche's method penalty parameter, larger values result in less erroneous penetration. Generally, should be ~100 times the Young's modulus.
  - `1.0`

* - `-[prefix]_friction_type [none | coulomb | threlfall]`
  - Specify the type of friction to use, see {ref}`friction-models` for details on available models.
  - `none`

* - `-[prefix]_friction_kinetic [non-negative real]`
  - Kinetic (sliding) coefficient of friction $\mu_k$, unitless.
  - `0.1`

* - `-[prefix]_friction_static [non-negative real]`
  - Static (slipping) coefficient of friction $\mu_s > \mu_k$, unitless. Not currently used.
  - `0.15`

* - `-[prefix]_friction_viscous [non-negative real]`
  - Viscous damping coefficient $F_v$, units of $\frac{\mathrm{N}\cdot\mathrm{s}}{\mathrm{m}}$.
  - `0.0`

* - `-[prefix]_friction_tolerance_velocity [positive real]`
  - Tolerance velocity $v_0$ at which full frictional slipping is permitted, units of $\frac{\mathrm{m}}{\mathrm{s}^2}$. Only used by Threlfall friction model.
  - `0.05`
:::

(using-contact-shapes)=

#### Shape-Specific Options
Each shape has different options for specifying its initial position and geometric properties.

(using-contact-shape-platen)=

##### Platen

(table_contact_platen_options)=

:::{list-table} Platen Shape Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-[prefix]_shape platen`
  - Needed to set the shape to a platen. For historical reasons, `platen` is the default shape.
  - `platen`

* - `-[prefix]_center [cx,cy,cz]`
  - Specify the center of the platen, with components given with respect to the global coordinate system
  -

* - `-[prefix]_normal [nx,ny,nz]`
  - Specify the exterior normal to the platen, with components given with respect to the global coordinate system.
    This vector should point toward the face `[contact_face]`, although it need not be parallel to the face normal.
  -

* - `-[prefix]_distance [dist][,dist2,dist3,...]`
  - Distance(s) of the half-plane along the specified normal vector.
    If more than one distance is specified, transition times must be provided as `-[prefix]_times [t1,t2,...]`.
    This is equivalent to specifying the translation `dist*n_x,dist*n_y,dist*n_z` with `-[prefix]_translate`.
    Only supported for platen contact BCs.
  - `0.0`
:::

##### Cylinder

(table_contact_cylinder_options)=

:::{list-table} Cylinder Shape Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-[prefix]_shape cylinder`
  - Needed to set the shape to a cylinder. Required.
  - `platen`

* - `-[prefix]_center [cx,cy,cz]`
  - Specify a point on the center line of the cylinder, with components given with respect to the global coordinate system. Required.
  -

* - `-[prefix]_axis [ax,ay,az]`
  - Specify the direction vector of the linear axis of the cylinder.
    This vector will be automatically normalized. Required.
  -

* - `-[prefix]_radius [float]`
  - Radius of the cylinder. Required.
  -
:::

:::{note}
  It is often easier to specify command-line arguments in YAML format, since repeated prefixes can be replaced with a YAML dictionary.
  For example, consider the following command-line arguments specifying a platen boundary condition named `top` with friction on face with attribute label value `1`:
  ```
  -bc_contact top -bc_contact_top_label_value 1 -bc_contact_top_shape platen -bc_contact_top_normal -1,0,0 -bc_contact_top_center 0,0,0 -bc_contact_top_distance 0.1 -bc_contact_top_gamma 500 -bc_contact_top_friction_type coulomb -bc_contact_top_friction_kinetic 0.4 -bc_contact_top_friction_viscous 0.01
  ```
  Alternatively, the same options specified via a YAML options file are far more readable:
  ```yaml
  # platen_options.yml
  bc:
    platen: top
    platen_top:
      label:
        value: 1
      shape: platen
      normal: -1,0,0
      center: 0,0,0
      distance: 0.1
      gamma: 500
      friction:
        type: coulomb
        kinetic: 0.4
        viscous: 0.01
  ```
  The options can then be passed to Ratel via the `-options_file platen_options.yml` command.
:::

## Miscellaneous Utilities
### Mesh Remapping for MPM
For MPM problems which are defined by uniaxial compression or tension enforced with clamp boundaries and slip boundary conditions on all other faces,
e.g. confined compression, it is convenient to remesh the background mesh at each timestep to prevent cells from becoming empty.

This remapping can be specified via the following options:

:::{list-table} Mesh Remapping Options
:header-rows: 1
:widths: 40 50 10
* - Option
  - Description
  - Default Value

* - `-remap_direction [x | y | z]`
  - The coordinate direction $d$ to rescale, provided as a character `x`, `y`, or `z`
  -

* - `-remap_scale [positive real]`
  - Scaling factor $s\in(0,\infty]$ along the provided direction.
  -
:::

The implemented functionality is limited to scaling linearly along a single Cartesian dimension.
More complicated remeshings are possible, but often require elastic smoothing of the mesh.
For details of how the remapping is defined, see the dropdown below.

:::{dropdown} Mesh Remapping Formulation
Let $d \in \{0,1,2\}$ be the direction and $z_{\min} = \min(x[d])$ ($z_{\max} = \max(x[d]$) be the minimum (maximum) initial node coordinate value in the given direction. Define the initial height $h_0 = z_{\max} - z_{\min}$ and base coordinate value $z_0 = z_{\min}$.
Given a scaling factor $s\in(0,\infty]$, define the final height by $h_f = s\cdot h_0$.

Assume that the height of the mesh can be parametrized in time by the linear interpolation:

$$
h(t) = t\cdot h_f + (1-t)\cdot h_0
$$

Assume that for a given node with coordinates $\bm x_k$, the $d$ coordinate $z_k$ can be parametrized in time by the function

$$
z_k(t) = z_0 + \alpha_k \cdot h(t)
$$

where $\alpha_k$ depends on the coordinate position and global height in the previous timestep,

$$
\alpha_k = \frac{z_k(t_{n-1}) - z_0}{h(t_{n-1}}.
$$

This characterization is convenient, as the new coordinate at time $t_n$ can be written in terms of the previous coordinate and height at time $t_{n-1}$, $z_k(t_{n-1})$ and $h(t_{n-1})$ and a set of time-independent constants: $d, z_0, h_0, h_f$.
:::

Ratel includes an example of this functionality in action:
:::{dropdown} `ex02-quasistatic-elasticity-mpm-neo-hookean-current-remap-verification.yml`

The following example performs uniaxial compression of the unit cube on the `x+` face.
Importantly, note that for this problem, compression of $0.25$ units corresponds to a scale factor of $s=0.75$.
See the emphasized lines below.

To run this problem in the Ratel directory, use the command-line options
```console
./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-mpm-neo-hookean-damage-current-remap-verification.yml -ts_monitor_swarm :swarm.xmf
```
The resulting output can be opened using Paraview with the command
```
paraview swarm.xmf
```


:::{literalinclude} /examples/ymls/ex02-quasistatic-elasticity-mpm-neo-hookean-damage-current-remap-verification.yml
:language: yaml
:emphasize-lines: 42-44,50-54

:::