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:

Table 2 Mandatory Runtime Options

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:

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

$ ./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 -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/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: 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

Table 3 Additional Runtime Options

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
-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
-bc_clamp [int list]	List of face sets on which to displace by -bc_clamp_[facenumber]_translate [x,y,z][,x2,y2,z2,...] and/or bc_clamp_[facenumber]_rotate [rx,ry,rz,c_0,c_1][,r2x,r2y,r2z,c2_0,c2_1,...]. If more than one translation and/or rotation is specified, transition times must be provided as bc_clamp_[facenumber]_times [t1,t2,...]. If t1 > 0, a node with zero translation and rotation will be added implicitly at time 0. Note: The default for a clamped face is zero displacement.
-bc_clamp_[facenumber]_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 bc_clamp_[facenumber]_times [t1,t2,...].	0,0,0
-bc_clamp_[facenumber]_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 bc_clamp_[facenumber]_times [t1,t2,...].	0,0,0,0,0
-bc_clamp_[facenumber]_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
-bc_clamp_[facenumber]_interp_type	Interpolation type between specified rotations/translations (none or linear).	linear
-bc_clamp_field_[fieldnumber] [int list]	List of face sets on which to displace by -bc_clamp_field_[fieldnumber]_face_[facenumber]_translate [x,y,z] and/or bc_clamp_field_[fieldnumber]_face_[facenumber]_rotate [rx,ry,rz,c_0,c_1] for a single solution field. Also supports lists of rotations or translations, see details in bc_clamp above. Note: The default for a clamped face is zero displacement.
-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_slip [int list]	List of face sets on which to set slip boundary conditions for the components -bc_slip_[facenumber]_components [int list]
-bc_slip_[facenumber]_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.
-bc_slip_[facenumber]_translate [c1,...,ck][c1_2,...,ck_2,...,c1_n,...,ck_n]	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 bc_slip_[facenumber]_times [t1,t2,...].	Zero vector of length \(k\)
-bc_slip_[facenumber]_times [t1,t2,...]	Transition times between each rigid 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
-bc_slip_[facenumber]_interp_type	Interpolation type between specified translations (none or linear).	linear
-bc_slip_field_[fieldnumber] [int list]	List of face sets on which to set slip boundary conditions for the components -bc_slip_field_[fieldnumber]_face_[facenumber]_components [int list] for a single solution field. Also supports lists of translations, see details in bc_slip above.
-bc_traction [int list]	List of face sets on which to set traction boundary conditions with the traction vector(s) -bc_traction_[facenumber] [tx,ty,tz][,t2x,t2y,t2z,...].
-bc_traction_[facenumber] [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 bc_traction_[facenumber]_times [t1,t2,...].
-bc_traction_[facenumber]_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
-bc_traction_[facenumber]_interp_type	Interpolation type between specified traction vectors (none or linear).	linear
-bc_platen [int list]	List of face sets on which to set platen (half-plane) contact boundary conditions. This boundary condition is performed using Nitsche’s method.
-bc_platen_[facenumber]_center [cx,cy,cz]	Specify the center of the platen, with components given with respect to the global coordinate system
-bc_platen_[facenumber]_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 [facenumber].
-bc_platen_[facenumber]_distance [dist][,dist2,dist3,...]	Distance(s) of the half-plane along the specified normal vector. In the context of timestepping, the speed of the platen depends on the value of -bc_platen_[facenumber]_interp_type. If more than one distance is specified, transition times must be provided as bc_platen_[facenumber]_times [t1,t2,...].	0.0
-bc_platen_[facenumber]_times [t1,t2,...]	Transition times between each platen distance value, primarily used to control the platen 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 distance is provided, the time value defaults to 1; i.e., the platen will move the prescribed distance in 1 time unit.	1
-bc_platen_[facenumber]_interp_type	Interpolation type between specified distance values (none or linear). Note: none results in piecewise-constant displacement while linear results in piecewise-linear displacement and piecewise-constant velocity.	linear
-bc_platen_[facenumber]_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
-bc_platen_[facenumber]_f [non-negative real]	Coefficient of friction, or 0 for frictionless contact	0
bc_pressure [int list]	List of face sets on which to set pressure boundary conditions.
bc_pressure_[facenumber] [real]	Specify the pressure.
-model	Material model to use (elasticity-linear, elasticity-neo-hookean-current, elasticity-mooney-rivlin-initial, etc.)	elasticity-linear
-forcing	Forcing term option (none, constant, or mms)	none
-forcing_vec [real list]	Forcing vector	0,-1,0
-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	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 list]	List of face sets on which to compute surface force
-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
-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
-ts_monitor_surface_force	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
-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_checkpoint [filename]	Output binary file with solution checkpoints for continuation. Note: Binary viewer and extension are always used.
-ts_monitor_checkpoint_interval [int]	Number of time steps between checkpoint output	1
-checkpoint_final_solution [filename]	Output binary file with final solution checkpoint. Note: Automatic with ts_monitor_checkpoint.
-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:

$ ./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.

Table 4 Material Property Options

Option	Description	Model
-E [real]	Young’s modulus, \(E > 0\)	Neo-Hookean
-nu [real]	Poisson’s ratio, \(\nu < 0.5\)	Neo-Hookean or Mooney-Rivlin
-nu_smoother [real]	Poisson’s ratio for multigrid smoothers, \(\nu < 0.5\)	Neo-Hookean or Mooney-Rivlin
-mu_1 [real]	Mooney-Rivlin material constant, \(\mu_1 > 0\),	Mooney-Rivlin
-mu_2 [real]	Mooney-Rivlin material constant, \(\mu_2 > 0\)	Mooney-Rivlin
-yield_stress [real]	Initial yield stress threshold, \(\sigma_0 > 0\)	Linear Plasticity
-hardening_A [real]	Isotropic hardening parameter, \(A > 0\)	Linear Plasticity

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/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:

Table 5 Multiple Material Runtime Options

Option	Description	Default value
-material [string list]	List of names to use as labels for each material.
-{material name}_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	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}_E [real]	Young’s modulus, \(E > 0\)
-{material name}_nu [real]	Poisson’s ratio, \(\nu < 0.5\)
-{material name}_nu_smoother [real]	Poisson’s ratio for multigrid smoothers, \(\nu < 0.5\)

An example of specifying a two material quasistatic solve with YAML is provided in 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 -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.

Table 6 (Non)dimensionalization options

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.

Table 7 Characteristic units for metals

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, \(\operatorname{trace} \bm{E}\), \(\operatorname{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.

Table 8 Diagnostic quantities

Quantity	Linear elasticity	Neo-Hookean hyperelasticity
Cauchy stress tensor	\(\bm{\sigma}\)	\(\bm{\sigma} = \bm{F} \bm{S} \bm{F}^T / J \)
Pressure	\(\lambda \operatorname{trace} \bm{\epsilon}\)	\(\frac{\lambda}{2J} \left(J^2 - 1\right) \)
Volumetric strain	\(\operatorname{trace} \bm{\epsilon}\)	\(\operatorname{trace} \bm{E}\)
\(\operatorname{trace} \bm{E}^2\)	\(\operatorname{trace} \bm{\epsilon}^2\)	\(\operatorname{trace} \bm{E}^2\)
\( J \)	\(1 + \operatorname{trace} \bm{\epsilon}\)	\( J \)
Strain energy density	\(\frac{\lambda}{2} (\operatorname{trace} \bm{\epsilon})^2 + \mu \bm{\epsilon} : \bm{\epsilon}\)	\(\frac{\lambda}{4}\left(J^2 - 1 - 2 \log J \right) + \mu \operatorname{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}}\)

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.