Using Ratel¶
The Ratel library includes support for processing inputs and handling commandline options. Most examples and applications using Ratel will inherit these options.
Command Line Options¶
Ratel is controlled via commandline 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:
Option 
Description 


Path to mesh file in any format supported by PETSc.
Alternatively, a builtin mesh, such as 

List of face sets on which to displace by 
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., downloadexodusii
for Exodus support.
Consider the specific example of the mesh seen below:
With the sidesets defined in the figure, we provide here an example of a minimal set of command line options:
$ ./bin/ratelquasistatic 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 dm_plex_filename
, the user may use a DMPlex box mesh by specifying dm_plex_box_faces [int list]
, dm_plex_box_upper [real list]
, and dm_plex_box_lower [real list]
.
As an alternative example exploiting dm_plex_box_faces
, we consider a 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:
$ ./bin/ratelquasistatic model elasticityneohookeaninitial 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: theta = (c_0 + c_1 * cx) * loadIncrement
where cx = kx * x + ky * y + kz * z
, with kx
, ky
, 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
Option 
Description 
Default value 


CEED resource specifier 


Filepath to yml file with runtime options for materials with one field 


Numerical method used to run simulations, either Finite Element Method ( 


Total number of material points, divided evenly between cells in the domain. Only used if 


Number of material points per cell, rounded up to the nearest cube. If not provided, 


Initial position of material points percell.
If 


Polynomial order of solution basis functions 


Polynomial order of solution basis functions for all fields 


Increased quadrature space order; final order given by 


Order for diagnostic values mesh 
Same value as multigrid fine level order specified via 

Geometry order for diagnostic values mesh 


Filepath to binary file holding restart information and vector 


List of face sets on which to displace by 


One or more vectors specifying rigid translation of the face.
Note: If more than one vector is specified, transition times must be provided as 


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 


Transition times between each rigid rotation/translation.
Note: If the first specified time is after 


Interpolation type between specified rotations/translations ( 


List of face sets on which to displace by 


List of face sets on which to displace by 


List of face sets on which to set Dirichlet boundary conditions to match a MMS solution 


List of face sets on which to set Dirichlet boundary conditions to match a MMS solution for a single solution field 


List of face sets on which to set Dirichlet boundary conditions to match a MMS solution for a single solution field 


List of face sets on which to set slip boundary conditions for the components 


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 


One or more vectors specifying rigid 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 
Zero vector of length \(k\) 

Transition times between each rigid translation.
Note: If the first specified time is after 


Interpolation type between specified translations ( 


List of face sets on which to set slip boundary conditions for the components 


List of face sets on which to set slip boundary conditions for the components 


List of face sets on which to set traction boundary conditions with the traction vector(s) 


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 


Transition times between each traction vector.
Note: If the first specified time is after 


Interpolation type between specified traction vectors ( 


List of labels for each platen (halfplane) contact boundary condition to apply. Examples: “top,bottom” or “1,2” This boundary condition is performed using Nitsche’s method. 


Solver method to use for platen boundary, either 


Value of 


Specify the center of the platen, with components given with respect to the global coordinate system 


Specify the exterior normal to the platen, with components given with respect to the global coordinate system.
This vector should point toward the face 


Distance(s) of the halfplane along the specified normal vector.
In the context of timestepping, the speed of the platen depends on the value of 


Transition times between each platen distance value, primarily used to control the platen velocity.
Note: If the first specified time is after 


Interpolation type between specified distance values ( 


Nitsche’s method penalty parameter, larger values result in less erroneous penetration. Generally, should be ~100 times the Young’s modulus. 


Coefficient of friction, or 


List of face sets on which to set pressure boundary conditions. 


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 


Transition times between each pressure value.
Note: If the first specified time is after 


Material model to use ( 


Forcing term option ( 


Userspecified acceleration vectors for applying body forces. Scaled by the mass computed with material density 


Transition times between each acceleration vector.
Note: If the first specified time is after 


Interpolation type between specified acceleration vectors ( 


Submatrix 


Submatrix 


PC type for linear solver 


PMultigrid coarsening to use ( 


Array of orders for each multigrid level, in ascending order; fine grid order specified by 


Polynomial order of coarse grid basis functions for materials with one field.
This is only used with 


List of named faces or face set label values on which to compute surface force and centroids. 


Face set label value for a given 


Optional bounding box specifying a subregion of a face set label value on which to compute surface forces and centroids. Useful for meshes which lack labeled faces. 


Increased quadrature space order for celltoface surface force computation; final order given by 


Expected strain energy, for testing 


Expected max displacement in a particular direction ( 


Expected surface force on face 


Expected surface forces on face 


Expected centroid of face 


Expected centroid of face 


Output final solution for viewing, ASCII format to STDOUT is used if no viewer is passed 


Output final MPM fields at material points for viewing, only 


Output computed strain energy on each time step 


Output final solution for viewing, ASCII format to STDOUT is used if no viewer is passed 


Output centroid displacements and accurate surface forces computed using the volumetric residual operator (reaction force) for faces given by 


Like 


Output centroid displacements and approximate surface forces computed using celltoface quadrature for faces given by 


Monitor swarm fields at each timestep, saving an XDMF file 


Output binary file with solution checkpoints for continuation. Note: Binary viewer and extension are always used. 


Output binary file with final solution checkpoint.
Note: Automatic with 


Number of time steps between monitor output, for any of 


View PETSc 


View PETSc 


View PETSc performance log 


Deprecated: Use 


View comprehensive information about runtime options 
To verify the convergence of the linear elasticity formulation on a given mesh with the method of manufactured solutions, run:
$ ./bin/ratelquasistatic dm_plex_filename [mesh] order [order] model elasticitylinear 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.
Option 
Description 
Model 


Young’s modulus, \(E > 0\) 
NeoHookean 

Poisson’s ratio, \(\nu < 0.5\) 
NeoHookean, MooneyRivlin or Ogden 

Poisson’s ratio for multigrid smoothers, \(\nu < 0.5\) 
NeoHookean, MooneyRivlin or Ogden 

MooneyRivlin material constant, \(\mu_1 > 0\), 
MooneyRivlin 

MooneyRivlin material constant, \(\mu_2 > 0\) 
MooneyRivlin 

Ogden 


Ogden material constant, \(2 \secondlame = \sum_{j=1}^N m_j \alpha_j,\) with \(m_j \alpha_j > 0\) 
Ogden 

Initial yield stress threshold, \(\sigma_0 > 0\) 
Linear Plasticity 

Isotropic hardening parameter, \(A > 0\) 
Linear Plasticity 

True (default): AT1 model, False: AT2 
Phase Field Modeling of Fracture 

True: hybrid formulation, False (default): nonhybrid (AT1 or AT2) model, 
Phase Field Modeling of Fracture 

True (default): full Jacobian, False: Jacobian diagonal terms only (for quasiNewton methods), 
Phase Field Modeling of Fracture 

Fracture_toughness expressed as critical energy release rate, \(fracture_{toughness} > 0\) 
Phase Field Modeling of Fracture 

Length scale parameter of phasefield model, \(characteristic_{length} > 0\) 
Phase Field Modeling of Fracture 

Residual stiffness in full fracture state (\(\phi = 1\)), \(residual_{stiffness} > 0\) 
Phase Field Modeling of Fracture 

Viscosity for viscous regularization of damage problem, \(damage_{viscosity} >= 0\) 
Phase Field Modeling of Fracture 
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 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 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 {material name}_
) to specify other aspects for each material at runtime.
We also specify, for each material, which domain label values to use with rod_label_value 3 binder_label_value 4
.
To define material parameters such as \(E\) and \(\nu\), we now use 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:
$ ./bin/ratelquasistatic material rod,binder rod_model elasticitymooneyrivlininitial rod_mu_1 0.5 rod_mu_2 0.5 rod_nu 0.4 binder_label_value 3 binder_model elasticityneohookeaninitial 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:
Option 
Description 
Default value 


List of names to use as labels for each material. 


Material to use ( 


Domain label specifying the type of volume to use for specifying materials. Optional. 


Domain value specifying the volume to use for a given material. 


Density of materal, by default in kg/m^3 


Young’s modulus, \(E > 0\) 


Poisson’s ratio, \(\nu < 0.5\) 


Poisson’s ratio for multigrid smoothers, \(\nu < 0.5\) 


Forcing term option ( 


Userspecified acceleration vectors for applying body forces.
Scaled by the mass computed with material density 


Transition times between each acceleration vector.
Note: If the first specified time is after 


Interpolation type between specified acceleration vectors ( 

An example of specifying a two material quasistatic solve with YAML is provided in examples/ex02quasistaticelasticitymultimaterial.yml
.
Algebraic Solvers¶
The examples are configured to use the following NewtonKrylovMultigrid method by default.
Newtontype 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; trymg_levels_ksp_max_it 2
, for example.) Application of the linear operators for all levels with order \(p > 1\) is performed matrixfree 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 commandline 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 commandline 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.
Option 
Description 
Default value 


1 meter in scaled length units 


1 second in scaled time units 


1 kilogram in scaled mass units 

For example, consider a problem involving metals subject to gravity.
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., 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, 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.
Quantity 
Linear elasticity 
NeoHookean 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\) 
\(J \rho\) 
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.