1# libCEED: Solid Mechanics Example 2 3This page provides a description of the solid mechanics example for the 4libCEED library, based on PETSc. 5PETSc v3.17 or a development version of PETSc at commit 0e95d842 or later is required. 6 7This code solves the steady-state static momentum balance equations using unstructured high-order finite/spectral element spatial discretizations. 8In this mini-app, we consider three formulations used in solid mechanics applications: linear elasticity, Neo-Hookean hyperelasticity at small strain, and Neo-Hookean hyperelasticity at finite strain. 9All three of these formulations are for compressible materials. 10 11Build by using: 12 13``` 14make 15``` 16 17and run with: 18 19``` 20./elasticity -mesh [.exo file] -degree [degree] -nu [nu] -E [E] [boundary options] -problem [problem type] -forcing [forcing] -ceed [ceed] 21``` 22 23## Runtime options 24 25% inclusion-solids-marker 26 27The elasticity mini-app is controlled via command-line options, the following of which are mandatory. 28 29:::{list-table} Mandatory Runtime Options 30:header-rows: 1 31:widths: 3 7 32 33* - Option 34 - Description 35* - `-mesh [filename]` 36 - Path to mesh file in any format supported by PETSc. 37* - `-degree [int]` 38 - Polynomial degree of the finite element basis 39* - `-E [real]` 40 - [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus), $E > 0$ 41* - `-nu [real]` 42 - [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu < 0.5$ 43* - `-bc_clamp [int list]` 44 - List of face sets on which to displace by `-bc_clamp_[facenumber]_translate [x,y,z]` 45 and/or `bc_clamp_[facenumber]_rotate [rx,ry,rz,c_0,c_1]`. Note: The default 46 for a clamped face is zero displacement. All displacement is with respect to 47 the initial configuration. 48* - `-bc_traction [int list]` 49 - List of face sets on which to set traction boundary conditions with the 50 traction vector `-bc_traction_[facenumber] [tx,ty,tz]` 51::: 52 53:::{note} 54This solver can use any mesh format that PETSc's `DMPlex` can read (Exodus, Gmsh, Med, etc.). 55Our 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]. 56Note that many mesh formats require PETSc to be configured appropriately; e.g., `--download-exodusii` for Exodus support. 57::: 58 59Consider the specific example of the mesh seen below: 60 61```{image} https://github.com/jeremylt/ceedSampleMeshes/raw/master/cylinderDiagram.png 62``` 63 64With the sidesets defined in the figure, we provide here an example of a minimal set of command line options: 65 66``` 67./elasticity -mesh [.exo file] -degree 4 -E 1e6 -nu 0.3 -bc_clamp 998,999 -bc_clamp_998_translate 0,-0.5,1 68``` 69 70In 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$. 71 72As an alternative to specifying a mesh with {code}`-mesh`, 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]`. 73 74As 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 all sides. Sides 1 through 6 are rotated around $x$-axis: 75 76``` 77./elasticity -problem FSInitial-NH1 -E 1 -nu 0.3 -num_steps 40 -snes_linesearch_type cp -dm_plex_box_faces 4,4,4 -bc_clamp 1,2,3,4,5,6 -bc_clamp_1_rotate 0,0,1,0,.3 -bc_clamp_2_rotate 0,0,1,0,.3 -bc_clamp_3_rotate 0,0,1,0,.3 -bc_clamp_4_rotate 0,0,1,0,.3 -bc_clamp_5_rotate 0,0,1,0,.3 -bc_clamp_6_rotate 0,0,1,0,.3 78``` 79 80:::{note} 81If 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. 82::: 83 84On 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. 85 86The command line options just shown are the minimum requirements to run the mini-app, but additional options may also be set as follows 87 88:::{list-table} Additional Runtime Options 89:header-rows: 1 90 91* - Option 92 - Description 93 - Default value 94 95* - `-ceed` 96 - CEED resource specifier 97 - `/cpu/self` 98 99* - `-qextra` 100 - Number of extra quadrature points 101 - `0` 102 103* - `-test` 104 - Run in test mode 105 - 106 107* - `-problem` 108 - Problem to solve (`Linear`, `SS-NH`, `FSInitial-NH1`, etc.) 109 - `Linear` 110 111* - `-forcing` 112 - Forcing term option (`none`, `constant`, or `mms`) 113 - `none` 114 115* - `-forcing_vec` 116 - Forcing vector 117 - `0,-1,0` 118 119* - `-multigrid` 120 - Multigrid coarsening to use (`logarithmic`, `uniform` or `none`) 121 - `logarithmic` 122 123* - `-nu_smoother [real]` 124 - Poisson's ratio for multigrid smoothers, $\nu < 0.5$ 125 - 126 127* - `-num_steps` 128 - Number of load increments for continuation method 129 - `1` if `Linear` else `10` 130 131* - `-view_soln` 132 - Output solution at each load increment for viewing 133 - 134 135* - `-view_final_soln` 136 - Output solution at final load increment for viewing 137 - 138 139* - `-snes_view` 140 - View PETSc `SNES` nonlinear solver configuration 141 - 142 143* - `-log_view` 144 - View PETSc performance log 145 - 146 147* - `-output_dir` 148 - Output directory 149 - `.` 150 151* - `-help` 152 - View comprehensive information about run-time options 153 - 154::: 155 156To verify the convergence of the linear elasticity formulation on a given mesh with the method of manufactured solutions, run: 157 158``` 159./elasticity -mesh [mesh] -degree [degree] -nu [nu] -E [E] -forcing mms 160``` 161 162This option attempts to recover a known solution from an analytically computed forcing term. 163 164### On algebraic solvers 165 166This mini-app is configured to use the following Newton-Krylov-Multigrid method by default. 167 168- Newton-type methods for the nonlinear solve, with the hyperelasticity models globalized using load increments. 169- Preconditioned conjugate gradients to solve the symmetric positive definite linear systems arising at each Newton step. 170- Preconditioning via $p$-version multigrid coarsening to linear elements, with algebraic multigrid (PETSc's `GAMG`) for the coarse solve. 171 The default smoother uses degree 3 Chebyshev with Jacobi preconditioning. 172 (Lower degree is often faster, albeit less robust; try {code}`-outer_mg_levels_ksp_max_it 2`, for example.) 173 Application of the linear operators for all levels with degree $p > 1$ is performed matrix-free using analytic Newton linearization, while the lowest order $p = 1$ operators are assembled explicitly (using coloring at present). 174 175Many related solvers can be implemented by composing PETSc command-line options. 176 177### Nondimensionalization 178 179Quantities such as the Young's modulus vary over many orders of magnitude, and thus can lead to poorly scaled equations. 180One can nondimensionalize the model by choosing an alternate system of units, such that displacements and residuals are of reasonable scales. 181 182:::{list-table} (Non)dimensionalization options 183:header-rows: 1 184 185* - Option 186 - Description 187 - Default value 188 189* - `-units_meter` 190 - 1 meter in scaled length units 191 - `1` 192 193* - `-units_second` 194 - 1 second in scaled time units 195 - `1` 196 197* - `-units_kilogram` 198 - 1 kilogram in scaled mass units 199 - `1` 200::: 201 202For example, consider a problem involving metals subject to gravity. 203 204:::{list-table} Characteristic units for metals 205:header-rows: 1 206 207* - Quantity 208 - Typical value in SI units 209 210* - Displacement, $\bm u$ 211 - $1 \,\mathrm{cm} = 10^{-2} \,\mathrm m$ 212 213* - Young's modulus, $E$ 214 - $10^{11} \,\mathrm{Pa} = 10^{11} \,\mathrm{kg}\, \mathrm{m}^{-1}\, \mathrm s^{-2}$ 215 216* - Body force (gravity) on volume, $\int \rho \bm g$ 217 - $5 \cdot 10^4 \,\mathrm{kg}\, \mathrm m^{-2} \, \mathrm s^{-2} \cdot (\text{volume} \, \mathrm m^3)$ 218::: 219 220One 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. 221This 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. 222 223### Diagnostic Quantities 224 225Diagnostic quantities for viewing are provided when the command line options for visualization output, {code}`-view_soln` or {code}`-view_final_soln` are used. 226The 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$, $\lvert J \rvert$, and strain energy density. 227The table below summarizes the formulations of each of these quantities for each problem type. 228 229:::{list-table} Diagnostic quantities 230 :header-rows: 1 231 232 * - Quantity 233 - Linear Elasticity 234 - Hyperelasticity, Small Strain 235 - Hyperelasticity, Finite Strain 236 237 * - Pressure 238 - $\lambda \operatorname{trace} \bm{\epsilon}$ 239 - $\lambda \log \operatorname{trace} \bm{\epsilon}$ 240 - $\lambda \log J$ 241 242 * - Volumetric Strain 243 - $\operatorname{trace} \bm{\epsilon}$ 244 - $\operatorname{trace} \bm{\epsilon}$ 245 - $\operatorname{trace} \bm{E}$ 246 247 * - $\operatorname{trace} \bm{E}^2$ 248 - $\operatorname{trace} \bm{\epsilon}^2$ 249 - $\operatorname{trace} \bm{\epsilon}^2$ 250 - $\operatorname{trace} \bm{E}^2$ 251 252 * - $\lvert J \rvert$ 253 - $1 + \operatorname{trace} \bm{\epsilon}$ 254 - $1 + \operatorname{trace} \bm{\epsilon}$ 255 - $\lvert J \rvert$ 256 257 * - Strain Energy Density 258 - $\frac{\lambda}{2} (\operatorname{trace} \bm{\epsilon})^2 + \mu \bm{\epsilon} : \bm{\epsilon}$ 259 - $\lambda (1 + \operatorname{trace} \bm{\epsilon}) (\log(1 + \operatorname{trace} \bm{\epsilon} ) - 1) + \mu \bm{\epsilon} : \bm{\epsilon}$ 260 - $\frac{\lambda}{2}(\log J)^2 + \mu \operatorname{trace} \bm{E} - \mu \log J$ 261::: 262 263[cubit]: https://cubit.sandia.gov/ 264[this repository]: https://github.com/jeremylt/ceedSampleMeshes 265
