- PorousFlowDictatorThe UserObject that holds the list of PorousFlow variable names
C++ Type:UserObjectName
Description:The UserObject that holds the list of PorousFlow variable names
 - drained_coefficientVolumetric coefficient of thermal expansion of the drained porous skeleton (ie the porous rock without fluid, or with a fluid that is free to move in and out of the rock)
C++ Type:double
Description:Volumetric coefficient of thermal expansion of the drained porous skeleton (ie the porous rock without fluid, or with a fluid that is free to move in and out of the rock)
 
PorousFlowConstantThermalExpansionCoefficient
Computes the effective thermal expansion coefficient, (biot_coeff - porosity) * drained_coefficient + porosity * fluid_coefficient.
This Material computes (1) where is the Biot Modulus, is the porosity, is the drained volumetric thermal expansion coefficient, and is the fluid thermal volumetric expansion coefficient.
This Material is designed to work with the PorousFlowFullySaturatedMassTimeDerivative Kernel.
This quantity is computed during the initial stage of the simulation and is kept fixed thereafter.
Input Parameters
- at_nodesFalseEvaluate Material properties at nodes instead of quadpoints
Default:False
C++ Type:bool
Options:
Description:Evaluate Material properties at nodes instead of quadpoints
 - biot_coefficient1Biot coefficient
Default:1
C++ Type:double
Options:
Description:Biot coefficient
 - blockThe list of block ids (SubdomainID) that this object will be applied
C++ Type:std::vector<SubdomainName>
Options:
Description:The list of block ids (SubdomainID) that this object will be applied
 - boundaryThe list of boundary IDs from the mesh where this boundary condition applies
C++ Type:std::vector<BoundaryName>
Options:
Description:The list of boundary IDs from the mesh where this boundary condition applies
 - computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Options:
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
 - constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
 - fluid_coefficient0.00021Volumetric coefficient of thermal expansion for the fluid
Default:0.00021
C++ Type:double
Options:
Description:Volumetric coefficient of thermal expansion for the fluid
 
Optional Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Options:
Description:Adds user-defined labels for accessing object parameters via control logic.
 - enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Options:
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
Options:
Description:Determines whether this object is calculated using an implicit or explicit form
 - seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Options:
Description:The seed for the master random number generator
 - use_displaced_meshFalseWhether 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:False
C++ Type:bool
Options:
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
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Options:
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
 - outputsnone Vector of output names were you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Options:
Description:Vector of output names were you would like to restrict the output of variables(s) associated with this object
 
Outputs Parameters
Input Files
- (modules/porous_flow/test/tests/actions/basicthm_th.i)
 - (modules/porous_flow/examples/tutorial/05_tabulated.i)
 - (modules/porous_flow/examples/tutorial/03.i)
 - (modules/porous_flow/examples/tutorial/04.i)
 - (modules/porous_flow/test/tests/actions/basicthm_thm.i)
 - (modules/porous_flow/test/tests/thm_rehbinder/fixed_outer.i)
 - (modules/porous_flow/test/tests/thm_rehbinder/fixed_outer_rz.i)
 - (modules/porous_flow/examples/tutorial/05.i)
 - (modules/porous_flow/test/tests/thm_rehbinder/free_outer.i)
 - (modules/porous_flow/test/tests/dirackernels/injection_production.i)
 - (modules/porous_flow/test/tests/jacobian/mass01_fully_saturated.i)
 
(modules/porous_flow/test/tests/actions/basicthm_th.i)
# PorousFlowBasicTHM action with coupling_type = ThermoHydroGenerator
# (no mechanical effects)
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 3
    xmax = 10
    ymax = 3
  []
  [aquifer]
    input = gen
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 1 0'
    top_right = '10 2 0'
  []
  [injection_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'injection_area'
    normal = '-1 0 0'
    input = 'aquifer'
  []
  [outflow_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'outflow_area'
    normal = '1 0 0'
    input = 'injection_area'
  []
  [rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caprock aquifer'
    input = 'outflow_area'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [porepressure]
    initial_condition = 1e6
  []
  [temperature]
    initial_condition = 293
    scaling = 1e-6
  []
[]
[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = simple_fluid
[]
[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1.5e6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 313
    boundary = injection_area
  []
  [constant_outflow_porepressure]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    boundary = outflow_area
    pt_vals = '0 1e9'
    multipliers = '0 1e9'
    flux_function = 1e-6
    PT_shift = 1e6
  []
  [constant_outflow_temperature]
    type = DirichletBC
    variable = temperature
    value = 293
    boundary = outflow_area
  []
[]
[Modules]
  [FluidProperties]
    [simple_fluid]
      type = SimpleFluidProperties
    []
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2e-7
    fluid_bulk_modulus = 1e7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1e-13 0 0   0 1e-13 0   0 0 1e-13'
  []
  [permeability_caprock]
    type = PorousFlowPermeabilityConst
    block = caprock
    permeability = '1e-15 0 0   0 1e-15 0   0 0 1e-15'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caprock aquifer'
  []
[]
[Preconditioning]
  [basic]
    type = SMP
    full = true
  []
[]
[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1e4
  dt = 1e3
  nl_abs_tol = 1e-15
  nl_rel_tol = 1e-14
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/05_tabulated.i)
# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [porepressure]
    initial_condition = 1E6
  []
  [temperature]
    initial_condition = 313
    scaling = 1E-8
  []
[]
[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = tabulated_water
[]
[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 2E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 333
    boundary = injection_area
  []
[]
[Modules]
  [FluidProperties]
    [true_water]
      type = Water97FluidProperties
    []
    [tabulated_water]
      type = TabulatedFluidProperties
      fp = true_water
      temperature_min = 275
      interpolated_properties = 'density viscosity enthalpy internal_energy'
      fluid_property_file = water97_tabulated.csv
    []
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
[]
[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 = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/03.i)
# Darcy flow with heat advection and conduction
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [porepressure]
  []
  [temperature]
    initial_condition = 293
    scaling = 1E-8
  []
[]
[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = the_simple_fluid
[]
[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 313
    boundary = injection_area
  []
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E9
      viscosity = 1.0E-3
      density0 = 1000.0
      thermal_expansion = 0.0002
      cp = 4194
      cv = 4186
      porepressure_coefficient = 0
    []
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
[]
[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 = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/04.i)
# Darcy flow with heat advection and conduction, and elasticity
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  []
[]
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
  biot_coefficient = 1.0
[]
[Variables]
  [porepressure]
  []
  [temperature]
    initial_condition = 293
    scaling = 1E-8
  []
  [disp_x]
    scaling = 1E-10
  []
  [disp_y]
    scaling = 1E-10
  []
  [disp_z]
    scaling = 1E-10
  []
[]
[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydroMechanical
  gravity = '0 0 0'
  fp = the_simple_fluid
  eigenstrain_names = thermal_contribution
  use_displaced_mesh = false
[]
[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 313
    boundary = injection_area
  []
  [roller_tmax]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = dmax
  []
  [roller_tmin]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = dmin
  []
  [roller_top_bottom]
    type = DirichletBC
    variable = disp_z
    value = 0
    boundary = 'top bottom'
  []
  [cavity_pressure_x]
    type = Pressure
    boundary = injection_area
    variable = disp_x
    component = 0
    factor = 1E6
    use_displaced_mesh = false
  []
  [cavity_pressure_y]
    type = Pressure
    boundary = injection_area
    variable = disp_y
    component = 1
    factor = 1E6
    use_displaced_mesh = false
  []
[]
[AuxVariables]
  [stress_rr]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_pp]
    family = MONOMIAL
    order = CONSTANT
  []
[]
[AuxKernels]
  [stress_rr]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_rr
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  []
  [stress_pp]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_pp
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  []
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E9
      viscosity = 1.0E-3
      density0 = 1000.0
      thermal_expansion = 0.0002
      cp = 4194
      cv = 4186
      porepressure_coefficient = 0
    []
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = thermal_contribution
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 0.001 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 293
  []
  [stress]
    type = ComputeLinearElasticStress
  []
[]
[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 = 1E6
  dt = 1E5
  nl_abs_tol = 1E-15
  nl_rel_tol = 1E-14
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/test/tests/actions/basicthm_thm.i)
# PorousFlowBasicTHM action with coupling_type = ThermoHydroMechanicalGenerator
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 3
    xmax = 10
    ymax = 3
  []
  [aquifer]
    input = gen
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 1 0'
    top_right = '10 2 0'
  []
  [injection_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'injection_area'
    normal = '-1 0 0'
    input = 'aquifer'
  []
  [outflow_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'outflow_area'
    normal = '1 0 0'
    input = 'injection_area'
  []
  [rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caprock aquifer'
    input = 'outflow_area'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
  displacements = 'disp_x disp_y'
  biot_coefficient = 1.0
[]
[Variables]
  [porepressure]
    initial_condition = 1e6
  []
  [temperature]
    initial_condition = 293
    scaling = 1e-6
  []
  [disp_x]
    scaling = 1e-6
  []
  [disp_y]
    scaling = 1e-6
  []
[]
[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydroMechanical
  gravity = '0 0 0'
  fp = simple_fluid
  eigenstrain_names = thermal_contribution
  use_displaced_mesh = false
  add_stress_aux = false
[]
[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1.5e6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 313
    boundary = injection_area
  []
  [constant_outflow_porepressure]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    boundary = outflow_area
    pt_vals = '0 1e9'
    multipliers = '0 1e9'
    flux_function = 1e-6
    PT_shift = 1e6
  []
  [constant_outflow_temperature]
    type = DirichletBC
    variable = temperature
    value = 293
    boundary = outflow_area
  []
  [top_bottom]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = 'top bottom'
  []
  [right]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = right
  []
[]
[Modules]
  [FluidProperties]
    [simple_fluid]
      type = SimpleFluidProperties
    []
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2e-7
    fluid_bulk_modulus = 1e7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1e-13 0 0   0 1e-13 0   0 0 1e-13'
  []
  [permeability_caprock]
    type = PorousFlowPermeabilityConst
    block = caprock
    permeability = '1e-15 0 0   0 1e-15 0   0 0 1e-15'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caprock aquifer'
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5e9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = thermal_contribution
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 0.001
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 293
  []
  [stress]
    type = ComputeLinearElasticStress
  []
[]
[Preconditioning]
  [basic]
    type = SMP
    full = true
  []
[]
[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1e4
  dt = 1e3
  nl_abs_tol = 1e-12
  nl_rel_tol = 1E-10
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/test/tests/thm_rehbinder/fixed_outer.i)
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 40
    nt = 16
    rmin = 0.1
    rmax = 1
    dmin = 0.0
    dmax = 90
    growth_r = 1.1
  []
  [make3D]
    input = annular
    type = MeshExtruderGenerator
    bottom_sideset = bottom
    top_sideset = top
    extrusion_vector = '0 0 1'
    num_layers = 1
  []
[]
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
  biot_coefficient = 1.0
[]
[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [porepressure]
  []
  [temperature]
  []
[]
[BCs]
  [plane_strain]
    type = DirichletBC
    variable = disp_z
    value = 0
    boundary = 'top bottom'
  []
  [ymin]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = dmin
  []
  [xmin]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = dmax
  []
  [cavity_temperature]
    type = DirichletBC
    variable = temperature
    value = 1000
    boundary = rmin
  []
  [cavity_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1E6
    boundary = rmin
  []
  [cavity_zero_effective_stress_x]
    type = Pressure
    component = 0
    variable = disp_x
    function = 1E6
    boundary = rmin
    use_displaced_mesh = false
  []
  [cavity_zero_effective_stress_y]
    type = Pressure
    component = 1
    variable = disp_y
    function = 1E6
    boundary = rmin
    use_displaced_mesh = false
  []
  [outer_temperature]
    type = DirichletBC
    variable = temperature
    value = 0
    boundary = rmax
  []
  [outer_pressure]
    type = DirichletBC
    variable = porepressure
    value = 0
    boundary = rmax
  []
  [fixed_outer_x]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = rmax
  []
  [fixed_outer_y]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = rmax
  []
[]
[AuxVariables]
  [stress_rr]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_pp]
    family = MONOMIAL
    order = CONSTANT
  []
[]
[AuxKernels]
  [stress_rr]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_rr
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  []
  [stress_pp]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_pp
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  []
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = SimpleFluidProperties
      thermal_expansion = 0.0
      bulk_modulus = 1E12
      viscosity = 1.0E-3
      density0 = 1000.0
      cv = 1000.0
      cp = 1000.0
      porepressure_coefficient = 0.0
    []
  []
[]
[PorousFlowBasicTHM]
  coupling_type = ThermoHydroMechanical
  multiply_by_density = false
  add_stress_aux = true
  porepressure = porepressure
  temperature = temperature
  eigenstrain_names = thermal_contribution
  gravity = '0 0 0'
  fp = the_simple_fluid
[]
[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1E10
    poissons_ratio = 0.2
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = thermal_contribution
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 1E-6
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 0.0
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [porosity]
    type = PorousFlowPorosityConst # only the initial value of this is ever used
    porosity = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    solid_bulk_compliance = 1E-10
    fluid_bulk_modulus = 1E12
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0   0 1E-12 0   0 0 1E-12'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    fluid_coefficient = 1E-6
    drained_coefficient = 1E-6
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '1E6 0 0  0 1E6 0  0 0 1E6'
  []
[]
[VectorPostprocessors]
  [P]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = porepressure
  []
  [T]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = temperature
  []
  [U]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = disp_x
  []
[]
[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
    petsc_options_value = 'gmres      asm      lu           1E-8'
  []
[]
[Executioner]
  type = Steady
  solve_type = Newton
[]
[Outputs]
  file_base = fixed_outer
  execute_on = timestep_end
  csv = true
[]
(modules/porous_flow/test/tests/thm_rehbinder/fixed_outer_rz.i)
# A version of fixed_outer.i that uses the RZ cylindrical coordinate system
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 40 # this is the r direction
  ny = 1 # this is the height direction
  xmin = 0.1
  xmax = 1
  bias_x = 1.1
  ymin = 0.0
  ymax = 1.0
[]
[Problem]
  coord_type = RZ
[]
[GlobalParams]
  displacements = 'disp_r disp_z'
  PorousFlowDictator = dictator
  biot_coefficient = 1.0
[]
[Variables]
  [disp_r]
  []
  [disp_z]
  []
  [porepressure]
  []
  [temperature]
  []
[]
[BCs]
  [plane_strain]
    type = DirichletBC
    variable = disp_z
    value = 0
    boundary = 'top bottom'
  []
  [cavity_temperature]
    type = DirichletBC
    variable = temperature
    value = 1000
    boundary = left
  []
  [cavity_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1E6
    boundary = left
  []
  [cavity_zero_effective_stress_x]
    type = Pressure
    component = 0
    variable = disp_r
    function = 1E6
    boundary = left
    use_displaced_mesh = false
  []
  [outer_temperature]
    type = DirichletBC
    variable = temperature
    value = 0
    boundary = right
  []
  [outer_pressure]
    type = DirichletBC
    variable = porepressure
    value = 0
    boundary = right
  []
  [fixed_outer_disp]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = right
  []
[]
[AuxVariables]
  [stress_rr]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_pp]
    family = MONOMIAL
    order = CONSTANT
  []
[]
[AuxKernels]
  [stress_rr]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_rr
    index_i = 0
    index_j = 0
  []
  [stress_pp] # hoop stress
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_pp
    index_i = 2
    index_j = 2
  []
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = SimpleFluidProperties
      thermal_expansion = 0.0
      bulk_modulus = 1E12
      viscosity = 1.0E-3
      density0 = 1000.0
      cv = 1000.0
      cp = 1000.0
      porepressure_coefficient = 0.0
    []
  []
[]
[PorousFlowBasicTHM]
  coupling_type = ThermoHydroMechanical
  multiply_by_density = false
  add_stress_aux = true
  porepressure = porepressure
  temperature = temperature
  eigenstrain_names = thermal_contribution
  gravity = '0 0 0'
  fp = the_simple_fluid
[]
[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1E10
    poissons_ratio = 0.2
  []
  [strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = thermal_contribution
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 1E-6
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 0.0
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [porosity]
    type = PorousFlowPorosityConst # only the initial value of this is ever used
    porosity = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    solid_bulk_compliance = 1E-10
    fluid_bulk_modulus = 1E12
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0   0 1E-12 0   0 0 1E-12' # note this is ordered: rr, zz, angle-angle
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    fluid_coefficient = 1E-6
    drained_coefficient = 1E-6
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '1E6 0 0  0 1E6 0  0 0 1E6' # note this is ordered: rr, zz, angle-angle
  []
[]
[VectorPostprocessors]
  [P]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = porepressure
  []
  [T]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = temperature
  []
  [U]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = disp_r
  []
[]
[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
    petsc_options_value = 'gmres      asm      lu           1E-8'
  []
[]
[Executioner]
  type = Steady
  solve_type = Newton
[]
[Outputs]
  file_base = fixed_outer_rz
  execute_on = timestep_end
  csv = true
[]
(modules/porous_flow/examples/tutorial/05.i)
# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [porepressure]
    initial_condition = 1E6
  []
  [temperature]
    initial_condition = 313
    scaling = 1E-8
  []
[]
[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = the_simple_fluid
[]
[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 2E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 333
    boundary = injection_area
  []
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = Water97FluidProperties
    []
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
[]
[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 = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/test/tests/thm_rehbinder/free_outer.i)
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 40
    nt = 16
    rmin = 0.1
    rmax = 1
    dmin = 0.0
    dmax = 90
    growth_r = 1.1
  []
  [make3D]
    input = annular
    type = MeshExtruderGenerator
    bottom_sideset = bottom
    top_sideset = top
    extrusion_vector = '0 0 1'
    num_layers = 1
  []
[]
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
  biot_coefficient = 1.0
[]
[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [porepressure]
  []
  [temperature]
  []
[]
[BCs]
  # sideset 1 = outer
  # sideset 2 = cavity
  # sideset 3 = ymin
  # sideset 4 = xmin
  [plane_strain]
    type = DirichletBC
    variable = disp_z
    value = 0
    boundary = 'top bottom'
  []
  [ymin]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = dmin
  []
  [xmin]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = dmax
  []
  [cavity_temperature]
    type = DirichletBC
    variable = temperature
    value = 1000
    boundary = rmin
  []
  [cavity_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1E6
    boundary = rmin
  []
  [cavity_zero_effective_stress_x]
    type = Pressure
    component = 0
    variable = disp_x
    function = 1E6
    boundary = rmin
    use_displaced_mesh = false
  []
  [cavity_zero_effective_stress_y]
    type = Pressure
    component = 1
    variable = disp_y
    function = 1E6
    boundary = rmin
    use_displaced_mesh = false
  []
  [outer_temperature]
    type = DirichletBC
    variable = temperature
    value = 0
    boundary = rmax
  []
  [outer_pressure]
    type = DirichletBC
    variable = porepressure
    value = 0
    boundary = rmax
  []
[]
[AuxVariables]
  [stress_rr]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_pp]
    family = MONOMIAL
    order = CONSTANT
  []
[]
[AuxKernels]
  [stress_rr]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_rr
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  []
  [stress_pp]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_pp
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  []
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = SimpleFluidProperties
      thermal_expansion = 0.0
      bulk_modulus = 1E12
      viscosity = 1.0E-3
      density0 = 1000.0
      cv = 1000.0
      cp = 1000.0
      porepressure_coefficient = 0.0
    []
  []
[]
[PorousFlowBasicTHM]
  coupling_type = ThermoHydroMechanical
  multiply_by_density = false
  add_stress_aux = true
  porepressure = porepressure
  temperature = temperature
  eigenstrain_names = thermal_contribution
  gravity = '0 0 0'
  fp = the_simple_fluid
[]
[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1E10
    poissons_ratio = 0.2
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = thermal_contribution
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 1E-6
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 0.0
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [porosity]
    type = PorousFlowPorosityConst # only the initial value of this is ever used
    porosity = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    solid_bulk_compliance = 1E-10
    fluid_bulk_modulus = 1E12
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0   0 1E-12 0   0 0 1E-12'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    fluid_coefficient = 1E-6
    drained_coefficient = 1E-6
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '1E6 0 0  0 1E6 0  0 0 1E6'
  []
[]
[VectorPostprocessors]
  [P]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = porepressure
  []
  [T]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = temperature
  []
  [U]
    type = LineValueSampler
    start_point = '0.1 0 0'
    end_point = '1.0 0 0'
    num_points = 10
    sort_by = x
    variable = disp_x
  []
[]
[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
    petsc_options_value = 'gmres      asm      lu           1E-8'
  []
[]
[Executioner]
  type = Steady
  solve_type = Newton
[]
[Outputs]
  file_base = free_outer
  execute_on = timestep_end
  csv = true
[]
(modules/porous_flow/test/tests/dirackernels/injection_production.i)
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 10
    ny = 10
    nz = 1
    xmin = -50
    xmax = 50
    ymin = -50
    ymax = 50
    zmin = 0
    zmax = 10
  []
  [central_nodes]
    input = gen
    type = ExtraNodesetGenerator
    new_boundary = central_nodes
    coord = '0 0 0; 0 0 10'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [porepressure]
    initial_condition = 20E6
  []
  [temperature]
    initial_condition = 400
    scaling = 1E-6 # fluid enthalpy is roughly 1E6
  []
[]
[BCs]
  [injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 300
    boundary = central_nodes
  []
[]
[DiracKernels]
  [fluid_injection]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = injected_mass
    point_file = injection.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 21E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = -1
  []
  [fluid_production]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = produced_mass
    point_file = production.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [remove_heat_at_production_well]
    type = PorousFlowPeacemanBorehole
    variable = temperature
    SumQuantityUO = produced_heat
    point_file = production.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    use_enthalpy = true
    character = 1
  []
[]
[UserObjects]
  [injected_mass]
    type = PorousFlowSumQuantity
  []
  [produced_mass]
    type = PorousFlowSumQuantity
  []
  [produced_heat]
    type = PorousFlowSumQuantity
  []
[]
[Postprocessors]
  [heat_joules_extracted_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_heat
  []
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = SimpleFluidProperties
      thermal_expansion = 2E-4
      bulk_modulus = 2E9
      viscosity = 1E-3
      density0 = 1000
      cv = 4000.0
      cp = 4000.0
    []
  []
[]
[PorousFlowUnsaturated]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = the_simple_fluid
[]
[Materials]
  [porosity]
    type = PorousFlowPorosityConst # only the initial value of this is ever used
    porosity = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    solid_bulk_compliance = 1E-10
    fluid_bulk_modulus = 2E9
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0   0 1E-12 0   0 0 1E-12'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    fluid_coefficient = 5E-6
    drained_coefficient = 2E-4
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '1 0 0  0 1 0  0 0 1'
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
[]
[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 = 2E6
  dt = 2E5
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/test/tests/jacobian/mass01_fully_saturated.i)
# FullySaturatedMassTimeDerivative
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
  displacements = 'disp_x disp_y disp_z'
[]
[Modules]
  [FluidProperties]
    [the_simple_fluid]
      type = SimpleFluidProperties
      thermal_expansion = 0.5
      bulk_modulus = 1.5
      density0 = 1.0
    []
  []
[]
[Variables]
  [pp]
  []
  [T]
  []
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]
[ICs]
  [disp_x]
    type = RandomIC
    min = -0.1
    max = 0.1
    variable = disp_x
  []
  [disp_y]
    type = RandomIC
    min = -0.1
    max = 0.1
    variable = disp_y
  []
  [disp_z]
    type = RandomIC
    min = -0.1
    max = 0.1
    variable = disp_z
  []
  [pp]
    type = RandomIC
    variable = pp
    min = 0
    max = 1
  []
  [T]
    type = RandomIC
    variable = T
    min = 0
    max = 1
  []
[]
[BCs]
  # necessary otherwise volumetric strain rate will be zero
  [disp_x]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = 'left right'
  []
  [disp_y]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = 'left right'
  []
  [disp_z]
    type = DirichletBC
    variable = disp_z
    value = 0
    boundary = 'left right'
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowFullySaturatedMassTimeDerivative
    variable = pp
    coupling_type = ThermoHydroMechanical
    biot_coefficient = 0.9
  []
  [dummyT]
    type = TimeDerivative
    variable = T
  []
  [grad_stress_x]
    type = StressDivergenceTensors
    variable = disp_x
    component = 0
  []
  [grad_stress_y]
    type = StressDivergenceTensors
    variable = disp_y
    component = 1
  []
  [grad_stress_z]
    type = StressDivergenceTensors
    variable = disp_z
    component = 2
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp disp_x disp_y disp_z T'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [simple1]
    type = TensorMechanicsPlasticSimpleTester
    a = 0
    b = 1
    strength = 1E20
    yield_function_tolerance = 1.0E-9
    internal_constraint_tolerance = 1.0E-9
  []
[]
[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    bulk_modulus = 2.0
    shear_modulus = 3.0
  []
  [strain]
    type = ComputeSmallStrain
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
  [temperature]
    type = PorousFlowTemperature
    temperature = T
  []
  [ppss]
    type = PorousFlow1PhaseFullySaturated
    porepressure = pp
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = the_simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst  # only the initial vaue of this is ever used
    porosity = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.9
    fluid_bulk_modulus = 1.5
    solid_bulk_compliance = 0.5
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.9
    fluid_coefficient = 0.5
    drained_coefficient = 0.4
  []
[]
[Preconditioning]
  [check]
    type = SMP
    full = true
    #petsc_options = '-snes_test_display'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
    petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
  []
[]
[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1
  end_time = 2
[]
[Outputs]
  exodus = false
[]