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