Stress Divergence RZ Tensors

Calculate stress divergence for an axisymmetric problem in cylindrical coordinates.

Description

The kernel StressDivergenceRZTensors solves the stress divergence equation for an Axisymmetric problem in the cylindrical coordinate system on a 2D mesh.

warning:Symmetry Assumed About the

var element = document.getElementById("moose-equation-2f443652-c7da-41fc-bb3f-36457e2dc6ab");katex.render("z", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});

-axis

The axis of symmetry must lie along the $var element = document.getElementById("moose-equation-64ab692c-8060-473a-b0f4-cb498739c2a1");katex.render("z", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ -axis in a $var element = document.getElementById("moose-equation-63f1a597-b6f5-4575-b78d-fc1a78f2df36");katex.render("\\left(r, z, \\theta \\right)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ or cylindrical coordinate system. This symmetry orientation is required for the calculation of the residual and of the jacobian, as defined in Eq. (1).

The StressDivergenceRZTensors kernel can be automatically created with the Solid Mechanics Physics. Use of the tensor mechanics quasi-static physics is recommended to ensure the consistent setting of the use_displaced_mesh parameter for the strain formulation selected. For a detailed explanation of the settings for _use_displaced_mesh_ in mechanics problems and the Solid Mechanics Physics usage, see the Introduction/Stress Divergence page.

Residual Calculation

The stress divergence kernel handles the calculation of the residual, $var element = document.getElementById("moose-equation-f251f963-efec-4738-8f66-49d056928f30");katex.render("\\mathbb{R}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , from the governing equation and the calculation of the Jacobian. From the strong form of the governing equation for mechanics, neglecting body forces, $var element = document.getElementById("moose-equation-408d0c12-57f0-4607-a31f-507cd78ec4a4");katex.render("\\nabla \\cdot \\sigma = 0", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ the weak form, using Galerkin's method and the Gauss divergence theorem, becomes $var element = document.getElementById("moose-equation-b7adf7ea-3bf6-4d37-bcd5-a8edd63f8c62");katex.render("\\int_S n \\cdot (\\sigma \\psi) dS - \\int_V \\sigma \\cdot \\nabla \\psi dV = 0", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ in which $var element = document.getElementById("moose-equation-a2cdab9f-cbfe-4041-86bd-29caac382c41");katex.render("\\psi", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the test function. The second term of the weak form equation is the residual contribution calculated by the stress divergence kernel.

The calculation of the Jacobian can be approximated with the elasticity tensor if the simulation solve type is JFNK:

$var element = document.getElementById("moose-equation-675ffb7f-7986-4761-ad46-699ea72df121");katex.render("\\frac{d(\\sigma_{ij} \\cdot d(\\psi) / dx_j)}{du_k} = \\frac{d(C_{ijmn} \\cdot du_m / dx_n \\cdot d\\psi /dx_j)}{du_k}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ which is nonzero for $var element = document.getElementById("moose-equation-c998f927-205a-4f26-8d91-9aaaff8beac5");katex.render("m == k", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .

If the solve type for the simulation is set to NEWTON the finite deformation Jacobian will need to be calculated. Set the parameter use_finite_deform_jacobian = true in this case.

note:Use of the Solid Mechanics QuasiStatic Physics Recommended

The use_displaced_mesh parameter must be set correcting to ensure consistency in the equilibrium equation: if the stress is calculated with respect to the deformed mesh, the test function gradients must also be calculated with respect to the deformed mesh. The Solid Mechanics QuasiStatic Physics is designed to automatically determine and set the parameter correctly for the selected strain formulation. We recommend that users employ the Solid Mechanics QuasiStatic Physics whenever possible to ensure consistency between the test function gradients and the strain formulation selected.

In cylindrical coordinates, the divergence of a rank-2 tensor includes mixed term contributions. In the axisymmetric model we assume symmetric loading conditions, in addition to the zero out-of-plane shear strains, so that the residual computation is simplified.

$(1)var element = document.getElementById("moose-equation-e33c130b-ebdb-4eeb-a40f-759575f64af5");katex.render(" \\begin{aligned} \\nabla \\sigma & = \\left[ \\frac{\\partial \\sigma_{rr}}{\\partial r} + \\frac{u_r}{X_r}\\sigma_{\\theta \\theta} + \\frac{\\partial \\sigma_{rz}}{\\partial z} \\right] \\hat{e}_r \\\\ & + \\left[ \\frac{\\partial \\sigma_{zz}}{\\partial z} + \\frac{\\partial \\sigma_{rz}}{\\partial r} \\right] \\hat{e}_z \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

The calculation of the Jacobian is similarly complex, requiring up to four terms in the calculation of the diagonal entries.

note:Notation Order Change

The axisymmetric system changes the order of the displacement vector from $var element = document.getElementById("moose-equation-2dd7ea95-164c-4ce4-8d1f-bb98a5b46e5d");katex.render("(u_r, u_{\\theta}, u_z)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , usually seen in textbooks, to $var element = document.getElementById("moose-equation-f0f35785-29a4-4eba-8ef7-816f0df90f33");katex.render("(u_r, u_z, u_{\\theta})", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . Take care to follow this convention in your input files and when adding extra stresses.

Example Input File

note:Use RZ Coordinate Type

The coordinate type in the Problem block of the input file must be set to COORD_TYPE = RZ.

Using the solid mechanics quasi-static physics, as shown

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [./all]
        strain = FINITE
        add_variables = true
        block = 1
      [../]
    []
  []
[]

the StressDivergenceRZTensors kernel will be automatically built when the coordinate system in the Problem block is specified for the axisymmetric RZ system,

[Problem]
  coord_type = RZ
[]

and only two displacement variables are provided:

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

Input Parameters

componentAn integer corresponding to the direction the variable this kernel acts in. (0 refers to the radial and 1 to the axial displacement.)
C++ Type:unsigned int
Controllable:No
Description:An integer corresponding to the direction the variable this kernel acts in. (0 refers to the radial and 1 to the axial displacement.)
displacementsThe string of displacements suitable for the problem statement
C++ Type:std::vector<VariableName>
Controllable:No
Description:The string of displacements suitable for the problem statement
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

base_nameMaterial property base name
C++ Type:std::string
Controllable:No
Description:Material property base name
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
coupled_variables
C++ Type:std::vector<VariableName>
Controllable:No
eigenstrain_namesList of eigenstrains used in the strain calculation. Used for computing their derivatives for off-diagonal Jacobian terms.
C++ Type:std::vector<MaterialPropertyName>
Controllable:No
Description:List of eigenstrains used in the strain calculation. Used for computing their derivatives for off-diagonal Jacobian terms.
out_of_plane_directionzThe direction of the out_of_plane_strain variable used in the WeakPlaneStress kernel.
Default:z
C++ Type:MooseEnum
Options:x, y, z
Controllable:No
Description:The direction of the out_of_plane_strain variable used in the WeakPlaneStress kernel.
out_of_plane_strainThe name of the out_of_plane_strain variable used in the WeakPlaneStress kernel.
C++ Type:std::vector<VariableName>
Controllable:No
Description:The name of the out_of_plane_strain variable used in the WeakPlaneStress kernel.
prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
temperatureThe name of the temperature variable used in the ComputeThermalExpansionEigenstrain. (Not required for simulations without temperature coupling.)
C++ Type:std::vector<VariableName>
Controllable:No
Description:The name of the temperature variable used in the ComputeThermalExpansionEigenstrain. (Not required for simulations without temperature coupling.)
use_finite_deform_jacobianFalseJacobian for corotational finite strain
Default:False
C++ Type:bool
Controllable:No
Description:Jacobian for corotational finite strain
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
volumetric_locking_correctionFalseSet to false to turn off volumetric locking correction
Default:False
C++ Type:bool
Controllable:No
Description:Set to false to turn off volumetric locking correction

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Controllable:No
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Controllable:No
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshTrueWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:True
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

The stress divergence family of automatic differentiation kernels

ADStressDivergenceTensors - Stress divergence kernel for Cartesian coordinates Problem/coord_type=XYZ (default)
ADStressDivergenceRZTensors - Stress divergence kernel for 2D cylindrical coordinates Problem/coord_type=RZ
ADStressDivergenceRSphericalTensors - Stress divergence kernel for 1D spherical coordinates Problem/coord_type=RSPHERICAL

Input Files

(modules/porous_flow/examples/tutorial/11_2D.i)
(modules/solid_mechanics/test/tests/ad_elastic/rz_finite_elastic-noad.i)
(modules/solid_mechanics/test/tests/ad_elastic/rz_small_elastic-noad.i)
(modules/porous_flow/examples/thm_example/2D.i)
(modules/porous_flow/examples/thm_example/2D_c.i)
(modules/porous_flow/test/tests/mass_conservation/mass12.i)
(modules/solid_mechanics/test/tests/ad_elastic/rz_incremental_small_elastic-noad.i)
(modules/solid_mechanics/test/tests/jacobian/thermal_coupling_rz.i)
(modules/porous_flow/test/tests/mass_conservation/mass13.i)

(modules/solid_mechanics/test/tests/2D_geometries/2D-RZ_finiteStrain_test.i)

# Considers the mechanics solution for a thick spherical shell that is uniformly
# pressurized on the inner and outer surfaces, using 2D axisymmetric geometry.
# This test uses the strain calculator ComputeAxisymmetricRZFiniteStrain,
# which is generated through the use of the SolidMechanics QuasiStatic Physics.
#
# From Roark (Formulas for Stress and Strain, McGraw-Hill, 1975), the radially-dependent
# circumferential stress in a uniformly pressurized thick spherical shell is given by:
#
# S(r) = [ Pi[ri^3(2r^3+ro^3)] - Po[ro^3(2r^3+ri^3)] ] / [2r^3(ro^3-ri^3)]
#
#   where:
#          Pi = inner pressure
#          Po = outer pressure
#          ri = inner radius
#          ro = outer radius
#
# The tests assume an inner and outer radii of 5 and 10, with internal and external
# pressures of 100000 and 200000 at t = 1.0, respectively. The resulting compressive
# tangential stress is largest at the inner wall and, from the above equation, has a
# value of -271429.
#
# RESULTS are below. Since stresses are average element values, values for the
# edge element and one-element-in are used to extrapolate the stress to the
# inner surface. The vesrion of the tests that are checked use the coarsest meshes.
#
#  Mesh    Radial elem   S(edge elem)  S(one elem in)  S(extrap to surf)
# 1D-SPH
# 2D-RZ        12 (x10)    -265004      -254665        -270174
#  3D          12 (6x6)    -261880      -252811        -266415
#
# 1D-SPH
# 2D-RZ        48 (x10)    -269853      -266710        -271425
#  3D          48 (10x10)  -268522      -265653        -269957
#
# The numerical solution converges to the analytical solution as the mesh is
# refined.

[Mesh]
  file = 2D-RZ_mesh.e
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Problem]
  coord_type = RZ
[]

[Physics/SolidMechanics/QuasiStatic]
  [./all]
    strain = FINITE
    add_variables = true
    block = 1
  [../]
[]

[AuxVariables]
  [./stress_theta]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_theta]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./stress_theta]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    variable = stress_theta
    execute_on = timestep_end
  [../]
  [./strain_theta]
    type = RankTwoAux
    rank_two_tensor = total_strain
    index_i = 2
    index_j = 2
    variable = strain_theta
    execute_on = timestep_end
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e10
    poissons_ratio = 0.345
    block = 1
  [../]

  [./_elastic_strain]
    type = ComputeFiniteStrainElasticStress
    block = 1
  [../]
[]

[BCs]
# pin particle along symmetry planes
  [./no_disp_r]
    type = DirichletBC
    variable = disp_r
    boundary = xzero
    value = 0.0
  [../]

  [./no_disp_z]
    type = DirichletBC
    variable = disp_z
    boundary = yzero
    value = 0.0
  [../]

# exterior and internal pressures
  [./exterior_pressure_r]
    type = Pressure
    variable = disp_r
    boundary = outer
    function = '200000*t'
  [../]

 [./exterior_pressure_z]
    type = Pressure
    variable = disp_z
    boundary = outer
    function = '200000*t'
  [../]

  [./interior_pressure_r]
    type = Pressure
    variable = disp_r
    boundary = inner
    function = '100000*t'
  [../]

  [./interior_pressure_z]
    type = Pressure
    variable = disp_z
    boundary = inner
    function = '100000*t'
  [../]
[]

[Debug]
    show_var_residual_norms = true
[]

[Executioner]
  type = Transient

  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
  petsc_options_value = '  201               hypre    boomeramg      10'

  line_search = 'none'

  #Preconditioned JFNK (default)
  solve_type = 'PJFNK'

  nl_rel_tol = 5e-9
  nl_abs_tol = 1e-10
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 50

  start_time = 0.0
  end_time = 0.2
  dt = 0.1
[]

[Postprocessors]
  [./strainTheta]
    type = ElementAverageValue
    variable = strain_theta
  [../]
  [./stressTheta]
    type = ElementAverageValue
    variable = stress_theta
  [../]
  [./stressTheta_pt]
    type = PointValue
    point = '5.0 0.0 0.0'
    #bottom inside edge for comparison to theory; use csv = true
    variable = stress_theta
  [../]
[]

[Outputs]
  exodus = true
[]

(modules/solid_mechanics/test/tests/2D_geometries/2D-RZ_finiteStrain_test.i)

# Considers the mechanics solution for a thick spherical shell that is uniformly
# pressurized on the inner and outer surfaces, using 2D axisymmetric geometry.
# This test uses the strain calculator ComputeAxisymmetricRZFiniteStrain,
# which is generated through the use of the SolidMechanics QuasiStatic Physics.
#
# From Roark (Formulas for Stress and Strain, McGraw-Hill, 1975), the radially-dependent
# circumferential stress in a uniformly pressurized thick spherical shell is given by:
#
# S(r) = [ Pi[ri^3(2r^3+ro^3)] - Po[ro^3(2r^3+ri^3)] ] / [2r^3(ro^3-ri^3)]
#
#   where:
#          Pi = inner pressure
#          Po = outer pressure
#          ri = inner radius
#          ro = outer radius
#
# The tests assume an inner and outer radii of 5 and 10, with internal and external
# pressures of 100000 and 200000 at t = 1.0, respectively. The resulting compressive
# tangential stress is largest at the inner wall and, from the above equation, has a
# value of -271429.
#
# RESULTS are below. Since stresses are average element values, values for the
# edge element and one-element-in are used to extrapolate the stress to the
# inner surface. The vesrion of the tests that are checked use the coarsest meshes.
#
#  Mesh    Radial elem   S(edge elem)  S(one elem in)  S(extrap to surf)
# 1D-SPH
# 2D-RZ        12 (x10)    -265004      -254665        -270174
#  3D          12 (6x6)    -261880      -252811        -266415
#
# 1D-SPH
# 2D-RZ        48 (x10)    -269853      -266710        -271425
#  3D          48 (10x10)  -268522      -265653        -269957
#
# The numerical solution converges to the analytical solution as the mesh is
# refined.

[Mesh]
  file = 2D-RZ_mesh.e
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Problem]
  coord_type = RZ
[]

[Physics/SolidMechanics/QuasiStatic]
  [./all]
    strain = FINITE
    add_variables = true
    block = 1
  [../]
[]

[AuxVariables]
  [./stress_theta]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_theta]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./stress_theta]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    variable = stress_theta
    execute_on = timestep_end
  [../]
  [./strain_theta]
    type = RankTwoAux
    rank_two_tensor = total_strain
    index_i = 2
    index_j = 2
    variable = strain_theta
    execute_on = timestep_end
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e10
    poissons_ratio = 0.345
    block = 1
  [../]

  [./_elastic_strain]
    type = ComputeFiniteStrainElasticStress
    block = 1
  [../]
[]

[BCs]
# pin particle along symmetry planes
  [./no_disp_r]
    type = DirichletBC
    variable = disp_r
    boundary = xzero
    value = 0.0
  [../]

  [./no_disp_z]
    type = DirichletBC
    variable = disp_z
    boundary = yzero
    value = 0.0
  [../]

# exterior and internal pressures
  [./exterior_pressure_r]
    type = Pressure
    variable = disp_r
    boundary = outer
    function = '200000*t'
  [../]

 [./exterior_pressure_z]
    type = Pressure
    variable = disp_z
    boundary = outer
    function = '200000*t'
  [../]

  [./interior_pressure_r]
    type = Pressure
    variable = disp_r
    boundary = inner
    function = '100000*t'
  [../]

  [./interior_pressure_z]
    type = Pressure
    variable = disp_z
    boundary = inner
    function = '100000*t'
  [../]
[]

[Debug]
    show_var_residual_norms = true
[]

[Executioner]
  type = Transient

  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
  petsc_options_value = '  201               hypre    boomeramg      10'

  line_search = 'none'

  #Preconditioned JFNK (default)
  solve_type = 'PJFNK'

  nl_rel_tol = 5e-9
  nl_abs_tol = 1e-10
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 50

  start_time = 0.0
  end_time = 0.2
  dt = 0.1
[]

[Postprocessors]
  [./strainTheta]
    type = ElementAverageValue
    variable = strain_theta
  [../]
  [./stressTheta]
    type = ElementAverageValue
    variable = stress_theta
  [../]
  [./stressTheta_pt]
    type = PointValue
    point = '5.0 0.0 0.0'
    #bottom inside edge for comparison to theory; use csv = true
    variable = stress_theta
  [../]
[]

[Outputs]
  exodus = true
[]

(modules/solid_mechanics/test/tests/2D_geometries/2D-RZ_finiteStrain_test.i)

# Considers the mechanics solution for a thick spherical shell that is uniformly
# pressurized on the inner and outer surfaces, using 2D axisymmetric geometry.
# This test uses the strain calculator ComputeAxisymmetricRZFiniteStrain,
# which is generated through the use of the SolidMechanics QuasiStatic Physics.
#
# From Roark (Formulas for Stress and Strain, McGraw-Hill, 1975), the radially-dependent
# circumferential stress in a uniformly pressurized thick spherical shell is given by:
#
# S(r) = [ Pi[ri^3(2r^3+ro^3)] - Po[ro^3(2r^3+ri^3)] ] / [2r^3(ro^3-ri^3)]
#
#   where:
#          Pi = inner pressure
#          Po = outer pressure
#          ri = inner radius
#          ro = outer radius
#
# The tests assume an inner and outer radii of 5 and 10, with internal and external
# pressures of 100000 and 200000 at t = 1.0, respectively. The resulting compressive
# tangential stress is largest at the inner wall and, from the above equation, has a
# value of -271429.
#
# RESULTS are below. Since stresses are average element values, values for the
# edge element and one-element-in are used to extrapolate the stress to the
# inner surface. The vesrion of the tests that are checked use the coarsest meshes.
#
#  Mesh    Radial elem   S(edge elem)  S(one elem in)  S(extrap to surf)
# 1D-SPH
# 2D-RZ        12 (x10)    -265004      -254665        -270174
#  3D          12 (6x6)    -261880      -252811        -266415
#
# 1D-SPH
# 2D-RZ        48 (x10)    -269853      -266710        -271425
#  3D          48 (10x10)  -268522      -265653        -269957
#
# The numerical solution converges to the analytical solution as the mesh is
# refined.

[Mesh]
  file = 2D-RZ_mesh.e
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Problem]
  coord_type = RZ
[]

[Physics/SolidMechanics/QuasiStatic]
  [./all]
    strain = FINITE
    add_variables = true
    block = 1
  [../]
[]

[AuxVariables]
  [./stress_theta]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./strain_theta]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./stress_theta]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    variable = stress_theta
    execute_on = timestep_end
  [../]
  [./strain_theta]
    type = RankTwoAux
    rank_two_tensor = total_strain
    index_i = 2
    index_j = 2
    variable = strain_theta
    execute_on = timestep_end
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e10
    poissons_ratio = 0.345
    block = 1
  [../]

  [./_elastic_strain]
    type = ComputeFiniteStrainElasticStress
    block = 1
  [../]
[]

[BCs]
# pin particle along symmetry planes
  [./no_disp_r]
    type = DirichletBC
    variable = disp_r
    boundary = xzero
    value = 0.0
  [../]

  [./no_disp_z]
    type = DirichletBC
    variable = disp_z
    boundary = yzero
    value = 0.0
  [../]

# exterior and internal pressures
  [./exterior_pressure_r]
    type = Pressure
    variable = disp_r
    boundary = outer
    function = '200000*t'
  [../]

 [./exterior_pressure_z]
    type = Pressure
    variable = disp_z
    boundary = outer
    function = '200000*t'
  [../]

  [./interior_pressure_r]
    type = Pressure
    variable = disp_r
    boundary = inner
    function = '100000*t'
  [../]

  [./interior_pressure_z]
    type = Pressure
    variable = disp_z
    boundary = inner
    function = '100000*t'
  [../]
[]

[Debug]
    show_var_residual_norms = true
[]

[Executioner]
  type = Transient

  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
  petsc_options_value = '  201               hypre    boomeramg      10'

  line_search = 'none'

  #Preconditioned JFNK (default)
  solve_type = 'PJFNK'

  nl_rel_tol = 5e-9
  nl_abs_tol = 1e-10
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 50

  start_time = 0.0
  end_time = 0.2
  dt = 0.1
[]

[Postprocessors]
  [./strainTheta]
    type = ElementAverageValue
    variable = strain_theta
  [../]
  [./stressTheta]
    type = ElementAverageValue
    variable = stress_theta
  [../]
  [./stressTheta_pt]
    type = PointValue
    point = '5.0 0.0 0.0'
    #bottom inside edge for comparison to theory; use csv = true
    variable = stress_theta
  [../]
[]

[Outputs]
  exodus = true
[]

(modules/porous_flow/examples/tutorial/11_2D.i)

# Two-phase borehole injection problem in RZ coordinates
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    xmin = 1.0
    xmax = 10
    bias_x = 1.4
    ny = 3
    ymin = -6
    ymax = 6
  []
  [aquifer]
    input = gen
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 -2 0'
    top_right = '10 2 0'
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x<1.0001'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
  coord_type = RZ
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_r'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  []
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
  gravity = '0 0 0'
  biot_coefficient = 1.0
  PorousFlowDictator = dictator
[]

[Variables]
  [pwater]
    initial_condition = 20E6
  []
  [pgas]
    initial_condition = 20.1E6
  []
  [T]
    initial_condition = 330
    scaling = 1E-5
  []
  [disp_r]
    scaling = 1E-5
  []
[]

[Kernels]
  [mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  []
  [vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    variable = pwater
  []
  [mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pgas
  []
  [flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  []
  [vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    variable = pgas
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = T
  []
  [advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  []
  [conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  []
  [vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    variable = T
  []
  [grad_stress_r]
    type = StressDivergenceRZTensors
    temperature = T
    variable = disp_r
    eigenstrain_names = thermal_contribution
    use_displaced_mesh = false
    component = 0
  []
  [poro_r]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_r
    use_displaced_mesh = false
    component = 0
  []
[]

[AuxVariables]
  [disp_z]
  []
  [effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  []
  [mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  []
  [mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  []
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
  [swater]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_rr]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_tt]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_zz]
    family = MONOMIAL
    order = CONSTANT
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [effective_fluid_pressure]
    type = ParsedAux
    coupled_variables = 'pwater pgas swater sgas'
    expression = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  []
  [swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [stress_rr_aux]
    type = RankTwoAux
    variable = stress_rr
    rank_two_tensor = stress
    index_i = 0
    index_j = 0
  []
  [stress_tt]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_tt
    index_i = 2
    index_j = 2
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 1
    index_j = 1
  []
  [porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = timestep_end
  []
[]

[BCs]
  [pinned_top_bottom_r]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = 'top bottom'
  []
  [cavity_pressure_r]
    type = Pressure
    boundary = injection_area
    variable = disp_r
    postprocessor = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh = false
  []

  [cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290 # injection temperature
    use_displaced_mesh = false
  []
  [constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  []

  [outer_water_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pwater
    fluid_phase = 0
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
  [outer_co2_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pgas
    fluid_phase = 1
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20.1E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedBicubicFluidProperties
    fp = true_water
    temperature_min = 275
    pressure_max = 1E8
    fluid_property_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties
  []
  [tabulated_co2]
    type = TabulatedBicubicFluidProperties
    fp = true_co2
    temperature_min = 275
    pressure_max = 1E8
    fluid_property_file = co2_tabulated_11.csv
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = T
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.1
    sum_s_res = 0.2
    phase = 0
  []
  [relperm_co2]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    thermal = true
    porosity_zero = 0.1
    reference_temperature = 330
    reference_porepressure = 20E6
    thermal_expansion_coeff = 15E-6 # volumetric
    solid_bulk = 8E9 # unimportant since biot = 1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityKozenyCarman
    block = aquifer
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-12
  []
  [permeability_caps]
    type = PorousFlowPermeabilityKozenyCarman
    block = caps
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-15
    k_anisotropy = '1 0 0  0 1 0  0 0 0.1'
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2 0 0  0 2 0  0 0 2'
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2300
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = 'thermal_contribution initial_stress'
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = T
    thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 330
  []
  [initial_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '20E6 0 0  0 20E6 0  0 0 20E6'
    eigenstrain_name = initial_stress
  []
  [stress]
    type = ComputeLinearElasticStress
  []

  [effective_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[Postprocessors]
  [effective_fluid_pressure_at_wellbore]
    type = PointValue
    variable = effective_fluid_pressure
    point = '1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
  [constrained_effective_fluid_pressure_at_wellbore]
    type = FunctionValuePostprocessor
    function = constrain_effective_fluid_pressure
    execute_on = timestep_begin
  []
[]

[Functions]
  [constrain_effective_fluid_pressure]
    type = ParsedFunction
    symbol_names = effective_fluid_pressure_at_wellbore
    symbol_values = effective_fluid_pressure_at_wellbore
    expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E3
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1E3
    growth_factor = 1.2
    optimal_iterations = 10
  []
  nl_abs_tol = 1E-7
[]

[Outputs]
  exodus = true
[]

(modules/solid_mechanics/test/tests/ad_elastic/rz_finite_elastic-noad.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 3
  ny = 3
[]

[Problem]
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Variables]
  # scale with one over Young's modulus
  [./disp_r]
    scaling = 1e-10
  [../]
  [./disp_z]
    scaling = 1e-10
  [../]
[]

[Kernels]
  [./stress_r]
    type = StressDivergenceRZTensors
    component = 0
    variable = disp_r
    use_displaced_mesh = true
  [../]
  [./stress_z]
    type = StressDivergenceRZTensors
    component = 1
    variable = disp_z
    use_displaced_mesh = true
  [../]
[]

[BCs]
  [./bottom]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0
  [../]
  [./axial]
    type = DirichletBC
    variable = disp_r
    boundary = left
    value = 0
  [../]
  [./rdisp]
    type = DirichletBC
    variable = disp_r
    boundary = right
    value = 0.1
  [../]
[]

[Materials]
  [./elasticity]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e10
  [../]
  [./strain]
    type = ComputeAxisymmetricRZFiniteStrain
  [../]
  [./stress]
    type = ComputeFiniteStrainElasticStress
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  dt = 0.05
  solve_type = 'NEWTON'

  petsc_options_iname = -pc_hypre_type
  petsc_options_value = boomeramg

  dtmin = 0.05
  num_steps = 1
[]

[Outputs]
  exodus = true
  file_base = rz_finite_elastic_out
[]

(modules/solid_mechanics/test/tests/ad_elastic/rz_small_elastic-noad.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 3
  ny = 3
[]

[Problem]
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Variables]
  # scale with one over Young's modulus
  [./disp_r]
    scaling = 1e-10
  [../]
  [./disp_z]
    scaling = 1e-10
  [../]
[]

[Kernels]
  [./stress_r]
    type = StressDivergenceRZTensors
    component = 0
    variable = disp_r
  [../]
  [./stress_z]
    type = StressDivergenceRZTensors
    component = 1
    variable = disp_z
  [../]
[]

[BCs]
  [./bottom]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0
  [../]
  [./axial]
    type = DirichletBC
    variable = disp_r
    boundary = left
    value = 0
  [../]
  [./rdisp]
    type = DirichletBC
    variable = disp_r
    boundary = right
    value = 0.1
  [../]
[]

[Materials]
  [./elasticity]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e10
  [../]
  [./strain]
    type = ComputeAxisymmetricRZSmallStrain
  [../]
  [./stress]
    type = ComputeLinearElasticStress
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  dt = 0.05
  solve_type = 'NEWTON'

  petsc_options_iname = -pc_hypre_type
  petsc_options_value = boomeramg

  dtmin = 0.05
  num_steps = 1
[]

[Outputs]
  exodus = true
  file_base = rz_small_elastic_out
[]

(modules/porous_flow/examples/thm_example/2D.i)

# Two phase, temperature-dependent, with mechanics, radial with fine mesh, constant injection of cold co2 into a overburden-reservoir-underburden containing mostly water
# species=0 is water
# species=1 is co2
# phase=0 is liquid, and since massfrac_ph0_sp0 = 1, this is all water
# phase=1 is gas, and since massfrac_ph1_sp0 = 0, this is all co2
#
# The mesh used below has very high resolution, so the simulation takes a long time to complete.
# Some suggested meshes of different resolution:
# nx=50, bias_x=1.2
# nx=100, bias_x=1.1
# nx=200, bias_x=1.05
# nx=400, bias_x=1.02
# nx=1000, bias_x=1.01
# nx=2000, bias_x=1.003
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2000
  bias_x = 1.003
  xmin = 0.1
  xmax = 5000
  ny = 1
  ymin = 0
  ymax = 11
[]

[Problem]
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
  PorousFlowDictator = dictator
  gravity = '0 0 0'
  biot_coefficient = 1.0
[]

[Variables]
  [pwater]
    initial_condition = 18.3e6
  []
  [sgas]
    initial_condition = 0.0
  []
  [temp]
    initial_condition = 358
  []
  [disp_r]
  []
[]

[AuxVariables]
  [rate]
  []
  [disp_z]
  []
  [massfrac_ph0_sp0]
    initial_condition = 1 # all H20 in phase=0
  []
  [massfrac_ph1_sp0]
    initial_condition = 0 # no H2O in phase=1
  []
  [pgas]
    family = MONOMIAL
    order = FIRST
  []
  [swater]
    family = MONOMIAL
    order = FIRST
  []
  [stress_rr]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_tt]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  []
  [mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = sgas
  []
  [flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = sgas
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = temp
  []
  [advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = temp
  []
  [conduction]
    type = PorousFlowExponentialDecay
    use_displaced_mesh = false
    variable = temp
    reference = 358
    rate = rate
  []
  [grad_stress_r]
    type = StressDivergenceRZTensors
    temperature = temp
    eigenstrain_names = thermal_contribution
    variable = disp_r
    use_displaced_mesh = false
    component = 0
  []
  [poro_r]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_r
    use_displaced_mesh = false
    component = 0
  []
[]

[AuxKernels]
  [rate]
    type = FunctionAux
    variable = rate
    execute_on = timestep_begin
    function = decay_rate
  []
  [pgas]
    type = PorousFlowPropertyAux
    property = pressure
    phase = 1
    variable = pgas
  []
  [swater]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 0
    variable = swater
  []
  [stress_rr]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_rr
    index_i = 0
    index_j = 0
  []
  [stress_tt]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_tt
    index_i = 2
    index_j = 2
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 1
    index_j = 1
  []
[]

[Functions]
  [decay_rate]
# Eqn(26) of the first paper of LaForce et al.
# Ka * (rho C)_a = 10056886.914
# h = 11
    type = ParsedFunction
    expression = 'sqrt(10056886.914/t)/11.0'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp pwater sgas disp_r'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
[]

[FluidProperties]
  [water]
    type = SimpleFluidProperties
    bulk_modulus = 2.27e14
    density0 = 970.0
    viscosity = 0.3394e-3
    cv = 4149.0
    cp = 4149.0
    porepressure_coefficient = 0.0
    thermal_expansion = 0
  []
  [co2]
    type = SimpleFluidProperties
    bulk_modulus = 2.27e14
    density0 = 516.48
    viscosity = 0.0393e-3
    cv = 2920.5
    cp = 2920.5
    porepressure_coefficient = 0.0
    thermal_expansion = 0
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow2PhasePS
    phase0_porepressure = pwater
    phase1_saturation = sgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = water
    phase = 0
  []
  [gas]
    type = PorousFlowSingleComponentFluid
    fp = co2
    phase = 1
  []
  [porosity_reservoir]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability_reservoir]
    type = PorousFlowPermeabilityConst
    permeability = '2e-12 0 0  0 0 0  0 0 0'
  []
  [relperm_liquid]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    phase = 0
    s_res = 0.200
    sum_s_res = 0.405
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityBC
    phase = 1
    s_res = 0.205
    sum_s_res = 0.405
    nw_phase = true
    lambda = 2
  []
  [thermal_conductivity_reservoir]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0 0 0  0 1.320 0  0 0 0'
    wet_thermal_conductivity = '0 0 0  0 3.083 0  0 0 0'
  []
  [internal_energy_reservoir]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2350.0
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    shear_modulus = 6.0E9
    poissons_ratio = 0.2
  []
  [strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = 'thermal_contribution ini_stress'
  []
  [ini_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '-12.8E6 0 0  0 -51.3E6 0  0 0 -12.8E6'
    eigenstrain_name = ini_stress
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temp
    stress_free_temperature = 358
    thermal_expansion_coeff = 5E-6
    eigenstrain_name = thermal_contribution
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[BCs]
  [outer_pressure_fixed]
    type = DirichletBC
    boundary = right
    value = 18.3e6
    variable = pwater
  []
  [outer_saturation_fixed]
    type = DirichletBC
    boundary = right
    value = 0.0
    variable = sgas
  []
  [outer_temp_fixed]
    type = DirichletBC
    boundary = right
    value = 358
    variable = temp
  []
  [fixed_outer_r]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = right
  []
  [co2_injection]
    type = PorousFlowSink
    boundary = left
    variable = sgas
    use_mobility = false
    use_relperm = false
    fluid_phase = 1
    flux_function = 'min(t/100.0,1)*(-2.294001475)' # 5.0E5 T/year = 15.855 kg/s, over area of 2Pi*0.1*11
  []
  [cold_co2]
    type = DirichletBC
    boundary = left
    variable = temp
    value = 294
  []
  [cavity_pressure_x]
    type = Pressure
    boundary = left
    variable = disp_r
    component = 0
    postprocessor = p_bh # note, this lags
    use_displaced_mesh = false
  []
[]

[Postprocessors]
  [p_bh]
    type = PointValue
    variable = pwater
    point = '0.1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
[]

[VectorPostprocessors]
  [ptsuss]
    type = LineValueSampler
    use_displaced_mesh = false
    start_point = '0.1 0 0'
    end_point = '5000 0 0'
    sort_by = x
    num_points = 50000
    outputs = csv
    variable = 'pwater temp sgas disp_r stress_rr stress_tt'
  []
[]

[Preconditioning]
  active = 'smp'
  [smp]
    type = SMP
    full = true
    #petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E2       1E-5        500'
  []
  [mumps]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package -pc_factor_shift_type -snes_rtol -snes_atol -snes_max_it'
    petsc_options_value = 'gmres      lu       mumps                         NONZERO               1E-5       1E2       50'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1.5768e8
  #dtmax = 1e6
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    growth_factor = 1.1
  []
[]

[Outputs]
  print_linear_residuals = false
  sync_times = '3600 86400 2.592E6 1.5768E8'
  perf_graph = true
  exodus = true
  [csv]
    type = CSV
    sync_only = true
  []
[]

(modules/porous_flow/examples/thm_example/2D_c.i)

# Two phase, temperature-dependent, with mechanics and chemistry, radial with fine mesh, constant injection of cold co2 into a overburden-reservoir-underburden containing mostly water
# species=0 is water
# species=1 is co2
# phase=0 is liquid, and since massfrac_ph0_sp0 = 1, this is all water
# phase=1 is gas, and since massfrac_ph1_sp0 = 0, this is all co2
#
# The mesh used below has very high resolution, so the simulation takes a long time to complete.
# Some suggested meshes of different resolution:
# nx=50, bias_x=1.2
# nx=100, bias_x=1.1
# nx=200, bias_x=1.05
# nx=400, bias_x=1.02
# nx=1000, bias_x=1.01
# nx=2000, bias_x=1.003
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2000
  bias_x = 1.003
  xmin = 0.1
  xmax = 5000
  ny = 1
  ymin = 0
  ymax = 11
[]

[Problem]
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
  PorousFlowDictator = dictator
  gravity = '0 0 0'
  biot_coefficient = 1.0
[]

[Variables]
  [pwater]
    initial_condition = 18.3e6
  []
  [sgas]
    initial_condition = 0.0
  []
  [temp]
    initial_condition = 358
  []
  [disp_r]
  []
[]

[AuxVariables]
  [rate]
  []
  [disp_z]
  []
  [massfrac_ph0_sp0]
    initial_condition = 1 # all H20 in phase=0
  []
  [massfrac_ph1_sp0]
    initial_condition = 0 # no H2O in phase=1
  []
  [pgas]
    family = MONOMIAL
    order = FIRST
  []
  [swater]
    family = MONOMIAL
    order = FIRST
  []
  [stress_rr]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_tt]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [mineral_conc_m3_per_m3]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = 0.1
  []
  [eqm_const]
    initial_condition = 0.0
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Kernels]
  [mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  []
  [mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = sgas
  []
  [flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = sgas
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = temp
  []
  [advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = temp
  []
  [conduction]
    type = PorousFlowExponentialDecay
    use_displaced_mesh = false
    variable = temp
    reference = 358
    rate = rate
  []
  [grad_stress_r]
    type = StressDivergenceRZTensors
    temperature = temp
    eigenstrain_names = thermal_contribution
    variable = disp_r
    use_displaced_mesh = false
    component = 0
  []
  [poro_r]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_r
    use_displaced_mesh = false
    component = 0
  []
[]

[AuxKernels]
  [rate]
    type = FunctionAux
    variable = rate
    execute_on = timestep_begin
    function = decay_rate
  []
  [pgas]
    type = PorousFlowPropertyAux
    property = pressure
    phase = 1
    variable = pgas
  []
  [swater]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 0
    variable = swater
  []
  [stress_rr]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_rr
    index_i = 0
    index_j = 0
  []
  [stress_tt]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_tt
    index_i = 2
    index_j = 2
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 1
    index_j = 1
  []
  [mineral]
    type = PorousFlowPropertyAux
    property = mineral_concentration
    mineral_species = 0
    variable = mineral_conc_m3_per_m3
  []
  [eqm_const_auxk]
    type = ParsedAux
    variable = eqm_const
    coupled_variables = temp
    expression = '(358 - temp) / (358 - 294)'
  []
  [porosity_auxk]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
  []
[]

[Functions]
  [decay_rate]
# Eqn(26) of the first paper of LaForce et al.
# Ka * (rho C)_a = 10056886.914
# h = 11
    type = ParsedFunction
    expression = 'sqrt(10056886.914/t)/11.0'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp pwater sgas disp_r'
    number_fluid_phases = 2
    number_fluid_components = 2
    number_aqueous_kinetic = 1
    aqueous_phase_number = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
[]

[FluidProperties]
  [water]
    type = SimpleFluidProperties
    bulk_modulus = 2.27e14
    density0 = 970.0
    viscosity = 0.3394e-3
    cv = 4149.0
    cp = 4149.0
    porepressure_coefficient = 0.0
    thermal_expansion = 0
  []
  [co2]
    type = SimpleFluidProperties
    bulk_modulus = 2.27e14
    density0 = 516.48
    viscosity = 0.0393e-3
    cv = 2920.5
    cp = 2920.5
    porepressure_coefficient = 0.0
    thermal_expansion = 0
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow2PhasePS
    phase0_porepressure = pwater
    phase1_saturation = sgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = water
    phase = 0
  []
  [gas]
    type = PorousFlowSingleComponentFluid
    fp = co2
    phase = 1
  []
  [porosity_reservoir]
    type = PorousFlowPorosity
    porosity_zero = 0.2
    chemical = true
    reference_chemistry = 0.1
    initial_mineral_concentrations = 0.1
  []
  [permeability_reservoir]
    type = PorousFlowPermeabilityConst
    permeability = '2e-12 0 0  0 0 0  0 0 0'
  []
  [relperm_liquid]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    phase = 0
    s_res = 0.200
    sum_s_res = 0.405
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityBC
    phase = 1
    s_res = 0.205
    sum_s_res = 0.405
    nw_phase = true
    lambda = 2
  []
  [thermal_conductivity_reservoir]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0 0 0  0 1.320 0  0 0 0'
    wet_thermal_conductivity = '0 0 0  0 3.083 0  0 0 0'
  []
  [internal_energy_reservoir]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2350.0
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    shear_modulus = 6.0E9
    poissons_ratio = 0.2
  []
  [strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = 'thermal_contribution ini_stress'
  []
  [ini_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '-12.8E6 0 0  0 -51.3E6 0  0 0 -12.8E6'
    eigenstrain_name = ini_stress
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temp
    stress_free_temperature = 358
    thermal_expansion_coeff = 5E-6
    eigenstrain_name = thermal_contribution
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
  [predis]
    type = PorousFlowAqueousPreDisChemistry
    num_reactions = 1
    primary_concentrations = 1.0 # fixed activity
    equilibrium_constants_as_log10 = true
    equilibrium_constants = eqm_const
    primary_activity_coefficients = 1.0 # fixed activity
    reactions = 1
    kinetic_rate_constant = 1E-6
    molar_volume = 1.0
    specific_reactive_surface_area = 1.0
    activation_energy = 0.0 # no Arrhenius
  []
  [mineral_conc]
    type = PorousFlowAqueousPreDisMineral
    initial_concentrations = 0.1
  []
  [predis_nodes]
    type = PorousFlowAqueousPreDisChemistry
    at_nodes = true
    num_reactions = 1
    primary_concentrations = 1.0 # fixed activity
    equilibrium_constants_as_log10 = true
    equilibrium_constants = eqm_const
    primary_activity_coefficients = 1.0 # fixed activity
    reactions = 1
    kinetic_rate_constant = 1E-6
    molar_volume = 1.0
    specific_reactive_surface_area = 1.0
    activation_energy = 0.0 # no Arrhenius
  []
  [mineral_conc_nodes]
    type = PorousFlowAqueousPreDisMineral
    at_nodes = true
    initial_concentrations = 0.1
  []
[]

[BCs]
  [outer_pressure_fixed]
    type = DirichletBC
    boundary = right
    value = 18.3e6
    variable = pwater
  []
  [outer_saturation_fixed]
    type = DirichletBC
    boundary = right
    value = 0.0
    variable = sgas
  []
  [outer_temp_fixed]
    type = DirichletBC
    boundary = right
    value = 358
    variable = temp
  []
  [fixed_outer_r]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = right
  []
  [co2_injection]
    type = PorousFlowSink
    boundary = left
    variable = sgas
    use_mobility = false
    use_relperm = false
    fluid_phase = 1
    flux_function = 'min(t/100.0,1)*(-2.294001475)' # 5.0E5 T/year = 15.855 kg/s, over area of 2Pi*0.1*11
  []
  [cold_co2]
    type = DirichletBC
    boundary = left
    variable = temp
    value = 294
  []
  [cavity_pressure_x]
    type = Pressure
    boundary = left
    variable = disp_r
    component = 0
    postprocessor = p_bh # note, this lags
    use_displaced_mesh = false
  []
[]

[Postprocessors]
  [p_bh]
    type = PointValue
    variable = pwater
    point = '0.1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
  [mineral_bh] # mineral concentration (m^3(mineral)/m^3(rock)) at the borehole
    type = PointValue
    variable = mineral_conc_m3_per_m3
    point = '0.1 0 0'
    use_displaced_mesh = false
  []
[]

[VectorPostprocessors]
  [ptsuss]
    type = LineValueSampler
    use_displaced_mesh = false
    start_point = '0.1 0 0'
    end_point = '5000 0 0'
    sort_by = x
    num_points = 50000
    outputs = csv
    variable = 'pwater temp sgas disp_r stress_rr stress_tt mineral_conc_m3_per_m3 porosity'
  []
[]

[Preconditioning]
  active = 'smp'
  [smp]
    type = SMP
    full = true
    #petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E2       1E-5        50'
  []
  [mumps]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package -pc_factor_shift_type -snes_rtol -snes_atol -snes_max_it'
    petsc_options_value = 'gmres      lu       mumps                         NONZERO               1E-5       1E2       50'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1.5768e8
  #dtmax = 1e6
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    growth_factor = 1.1
  []
[]

[Outputs]
  print_linear_residuals = false
  sync_times = '3600 86400 2.592E6 1.5768E8'
  perf_graph = true
  exodus = true
  [csv]
    type = CSV
    sync_only = true
  []
[]

(modules/porous_flow/test/tests/mass_conservation/mass12.i)

# The sample is an annulus in RZ coordinates.
# Roller BCs are applied to the r_min, r_max and bottom boundaries
# A constant displacement is applied to the top: disp_z = -0.01*t.
# There is no fluid flow.
# Fluid mass conservation is checked.
#
# Under these conditions
# fluid_mass = volume0 * rho0 * exp(P0/bulk) = pi*3 * 1 * exp(0.1/0.5) = 11.51145
# volume0 * rho0 * exp(P0/bulk) = volume * rho0 * exp(P/bulk), so
# P - P0 = bulk * log(volume0 / volume) = 0.5 * log(1 / (1 - 0.01*t))

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 1.0
  xmax = 2.0
  ymin = -0.5
  ymax = 0.5
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
  PorousFlowDictator = dictator
  block = 0
[]

[Variables]
  [disp_r]
  []
  [disp_z]
  []
  [porepressure]
    initial_condition = 0.1
  []
[]

[BCs]
  [bottom_roller]
    type = DirichletBC
    variable = disp_z
    value = 0
    boundary = bottom
  []
  [side_rollers]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = 'left right'
  []
  [top_move]
    type = FunctionDirichletBC
    variable = disp_z
    function = -0.01*t
    boundary = top
  []
[]

[Kernels]
  [grad_stress_r]
    type = StressDivergenceRZTensors
    variable = disp_r
    component = 0
  []
  [grad_stress_z]
    type = StressDivergenceRZTensors
    variable = disp_z
    component = 1
  []
  [poro_r]
    type = PorousFlowEffectiveStressCoupling
    biot_coefficient = 0.3
    variable = disp_r
    component = 0
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    biot_coefficient = 0.3
    variable = disp_z
    component = 1
  []
  [poro_vol_exp]
    type = PorousFlowMassVolumetricExpansion
    variable = porepressure
    fluid_component = 0
  []
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = porepressure
  []
[]

[AuxVariables]
  [stress_rr]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_rz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_tt]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [stress_rr]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_rr
    index_i = 0
    index_j = 0
  []
  [stress_rz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_rz
    index_i = 0
    index_j = 1
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 1
    index_j = 1
  []
  [stress_tt]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_tt
    index_i = 2
    index_j = 2
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 0.5
    density0 = 1
    thermal_expansion = 0
    viscosity = 1
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '1 1.5'
    # bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
    fill_method = symmetric_isotropic
  []
  [strain]
    type = ComputeAxisymmetricRZSmallStrain
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = porepressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.5 0 0   0 0.5 0   0 0 0.5'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_r disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
[]

[Postprocessors]
  [p0]
    type = PointValue
    outputs = 'console csv'
    execute_on = 'initial timestep_end'
    point = '1 0 0'
    variable = porepressure
  []
  [rdisp]
    type = PointValue
    outputs = csv
    point = '2 0 0'
    use_displaced_mesh = false
    variable = disp_r
  []
  [stress_rr]
    type = PointValue
    outputs = csv
    point = '1 0 0'
    variable = stress_rr
  []
  [stress_zz]
    type = PointValue
    outputs = csv
    point = '1 0 0'
    variable = stress_zz
  []
  [stress_tt]
    type = PointValue
    outputs = csv
    point = '1 0 0'
    variable = stress_tt
  []
  [fluid_mass]
    type = PorousFlowFluidMass
    fluid_component = 0
    execute_on = 'initial timestep_end'
    outputs = 'console csv'
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-14 1E-8 10000'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  start_time = 0
  end_time = 10
  dt = 2
[]

[Outputs]
  execute_on = 'initial timestep_end'
  [csv]
    type = CSV
  []
[]

(modules/solid_mechanics/test/tests/ad_elastic/rz_incremental_small_elastic-noad.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 3
  ny = 3
[]

[Problem]
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Variables]
  # scale with one over Young's modulus
  [./disp_r]
    scaling = 1e-10
  [../]
  [./disp_z]
    scaling = 1e-10
  [../]
[]

[Kernels]
  [./stress_r]
    type = StressDivergenceRZTensors
    component = 0
    variable = disp_r
  [../]
  [./stress_z]
    type = StressDivergenceRZTensors
    component = 1
    variable = disp_z
  [../]
[]

[BCs]
  [./bottom]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0
  [../]
  [./axial]
    type = DirichletBC
    variable = disp_r
    boundary = left
    value = 0
  [../]
  [./rdisp]
    type = DirichletBC
    variable = disp_r
    boundary = right
    value = 0.1
  [../]
[]

[Materials]
  [./elasticity]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.3
    youngs_modulus = 1e10
  [../]
  [./strain]
    type = ComputeAxisymmetricRZIncrementalStrain
  [../]
  [./stress]
    type = ComputeFiniteStrainElasticStress
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  dt = 0.05
  solve_type = 'NEWTON'

  petsc_options_iname = -pc_hypre_type
  petsc_options_value = boomeramg

  dtmin = 0.05
  num_steps = 1
[]

[Outputs]
  exodus = true
  file_base = rz_incremental_small_elastic_out
[]

(modules/solid_mechanics/test/tests/jacobian/thermal_coupling_rz.i)

# Thermal eigenstrain coupling
[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[Problem]
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Variables]
  [./disp_r]
  [../]
  [./disp_z]
  [../]
  [./temperature]
  [../]
[]

[Kernels]
  [./cx_elastic]
    type = StressDivergenceRZTensors
    variable = disp_r
    temperature = temperature
    eigenstrain_names = thermal_contribution
    use_displaced_mesh = false
    component = 0
  [../]
  [./cz_elastic]
    type = StressDivergenceRZTensors
    variable = disp_z
    temperature = temperature
    eigenstrain_names = thermal_contribution
    use_displaced_mesh = false
    component = 1
  [../]
  [./temperature]
    type = Diffusion
    variable = temperature
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 10.0
    poissons_ratio = 0.25
  [../]
  [./strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = thermal_contribution
  [../]
  [./thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 1.0E2
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 0.0
  [../]
  [./admissible]
    type = ComputeLinearElasticStress
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-snes_type'
    petsc_options_value = 'test'
  [../]
[]

[Executioner]
  solve_type = NEWTON
  end_time = 1
  dt = 1
  type = Transient
[]

(modules/porous_flow/test/tests/mass_conservation/mass13.i)

# The sample is an annulus in RZ coordinates.
# Roller BCs are applied to the rmin, top and bottom boundaries
# A constant displacement is applied to the outer boundary: disp_r = -0.01 * t * (r - rmin)/(rmax - rmin).
# There is no fluid flow.
# Fluid mass conservation is checked.
#
# The flag volumetric_locking_correction = true is set for the strain calculator,
# which ensures that the volumetric strain is uniform throughout the element
#
# Theoretically,
# volumetric_strain = volume / volume0 - 1 = ((rmax - 0.01*t)^2 - rmin^2) / (rmax^2 - rmin^2) - 1
# However, with ComputeAxisymmetricRZSmallStrain, strain_rr = -0.01 * t / (rmax - rmin)
# and strain_tt = disp_r / r = -0.01 * t * (1 - rmin / r_qp) / (rmax - rmin), where r_qp is the radius of the quadpoint
# With volumetric_locking_correction = true, r_qp = (rmax - rmin) / 2.
# The volumetric strain is
# epv = -0.01 * t * (2 - rmin / r_qp) / (rmax - rmin)
# and volume = volume0 * (1 + epv)
#
# Fluid conservation reads
# volume0 * rho0 * exp(P0/bulk) = volume * rho0 * exp(P/bulk), so
# P - P0 = bulk * log(volume0 / volume) = 0.5 * log(1 / (1 + epv))

# With rmax = 2 and rmin = 1
# fluid_mass = volume0 * rho0 * exp(P0/bulk) = pi*3 * 1 * exp(0.1/0.5) = 11.51145

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 1
  xmax = 2
  ymin = -0.5
  ymax = 0.5
  coord_type = RZ
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
  PorousFlowDictator = dictator
  block = 0
  biot_coefficient = 0.3
[]

[Variables]
  [disp_r]
  []
  [disp_z]
  []
  [porepressure]
    initial_condition = 0.1
  []
[]

[BCs]
  [plane_strain]
    type = DirichletBC
    variable = disp_z
    value = 0
    boundary = 'bottom top'
  []
  [rmin_fixed]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = left
  []
  [contract]
    type = FunctionDirichletBC
    variable = disp_r
    function = -0.01*t
    boundary = right
  []
[]

[Kernels]
  [grad_stress_r]
    type = StressDivergenceRZTensors
    variable = disp_r
    component = 0
  []
  [grad_stress_z]
    type = StressDivergenceRZTensors
    variable = disp_z
    component = 1
  []
  [poro_r]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_r
    component = 0
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_z
    component = 1
  []
  [poro_vol_exp]
    type = PorousFlowMassVolumetricExpansion
    variable = porepressure
    fluid_component = 0
  []
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = porepressure
  []
[]

[AuxVariables]
  [stress_rr]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_rz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_tt]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_rr]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_zz]
    order = CONSTANT
    family = MONOMIAL
  []
  [strain_tt]
    order = CONSTANT
    family = MONOMIAL
  []
  [vol_strain]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [stress_rr]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_rr
    index_i = 0
    index_j = 0
  []
  [stress_rz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_rz
    index_i = 0
    index_j = 1
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 1
    index_j = 1
  []
  [stress_tt]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_tt
    index_i = 2
    index_j = 2
  []
  [strain_rr]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_rr
    index_i = 0
    index_j = 0
  []
  [strain_zz]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_zz
    index_i = 1
    index_j = 1
  []
  [strain_tt]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = strain_tt
    index_i = 2
    index_j = 2
  []
  [vol_strain]
    type = MaterialRealAux
    property = PorousFlow_total_volumetric_strain_qp
    variable = vol_strain
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 0.5
    density0 = 1
    thermal_expansion = 0
    viscosity = 1
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '1 1.5'
    # bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
    fill_method = symmetric_isotropic
  []
  [strain]
    type = ComputeAxisymmetricRZSmallStrain
    volumetric_locking_correction = true # the strain will be the same at every qp of the element
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = porepressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.5 0 0   0 0.5 0   0 0 0.5'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_r disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
[]

[Postprocessors]
  [p0]
    type = PointValue
    outputs = 'console csv'
    execute_on = 'initial timestep_end'
    point = '1.0 0 0'
    variable = porepressure
  []
  [vol_strain]
    type = PointValue
    outputs = 'console csv'
    execute_on = 'initial timestep_end'
    point = '1 0 0'
    variable = vol_strain
  []
  [strain_rr]
    type = PointValue
    outputs = 'console csv'
    execute_on = 'initial timestep_end'
    point = '1 0 0'
    variable = strain_rr
  []
  [strain_zz]
    type = PointValue
    outputs = 'console csv'
    execute_on = 'initial timestep_end'
    point = '1 0 0'
    variable = strain_zz
  []
  [strain_tt]
    type = PointValue
    outputs = 'console csv'
    execute_on = 'initial timestep_end'
    point = '1 0 0'
    variable = strain_tt
  []
  [rdisp]
    type = PointValue
    outputs = csv
    point = '2 0 0'
    use_displaced_mesh = false
    variable = disp_r
  []
  [stress_rr]
    type = PointValue
    outputs = csv
    point = '1 0 0'
    variable = stress_rr
  []
  [stress_zz]
    type = PointValue
    outputs = csv
    point = '1 0 0'
    variable = stress_zz
  []
  [stress_tt]
    type = PointValue
    outputs = csv
    point = '1 0 0'
    variable = stress_tt
  []
  [fluid_mass]
    type = PorousFlowFluidMass
    fluid_component = 0
    execute_on = 'initial timestep_end'
    outputs = 'console csv'
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-14 1E-8 10000'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  start_time = 0
  end_time = 10
  dt = 2
[]

[Outputs]
  execute_on = 'initial timestep_end'
  [csv]
    type = CSV
  []
[]

Description
Residual Calculation
Example Input File
Input Parameters
The stress divergence family of automatic differentiation kernels
Input Files