PorousFlowPorosityLinear

This Material computes porosity (at the nodes or quadpoints, depending on the at_nodes flag) according to

$(1)var element = document.getElementById("moose-equation-8620f413-6427-44d3-a6b4-f3509c64f167");katex.render(" \\phi = \\phi_{\\mathrm{ref}} + A(P_{f} - P_{f\\ \\mathrm{ref}}) + B(T - T_{\\mathrm{ref}}) + C(\\epsilon^{\\mathrm{total}}_{ii} - \\epsilon^{\\mathrm{total}}_{ii\\ \\mathrm{ref}}) \\ .", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

Here:

$var element = document.getElementById("moose-equation-f2c7c5f8-cf9a-44b5-bf06-c17761bfd6a0");katex.render("\\phi_{\\mathrm{ref}}", 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 reference porosity, which is a constant-monomial variable or a real number
$var element = document.getElementById("moose-equation-76aa2a96-e888-4ce6-8b1a-9d562d0e1756");katex.render("A", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , $var element = document.getElementById("moose-equation-87bfd0ca-680f-48b0-9cd3-3c7e4789c41e");katex.render("B", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and $var element = document.getElementById("moose-equation-589002f4-b3b6-4e4a-8c53-9f1fa4276c97");katex.render("C", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ are real-valued coefficients
$var element = document.getElementById("moose-equation-2b7dc231-7ada-43bb-a361-0e9c5faed335");katex.render("P_{f}", 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 effective fluid pressure
$var element = document.getElementById("moose-equation-15724341-1345-4d20-9442-abf953f79f93");katex.render("P_{f\\ \\mathrm{ref}}", 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 reference effective fluid pressure, which is a coupled variable (usually linear-Lagrange) or a real number
$var element = document.getElementById("moose-equation-38093f93-034a-4765-9087-dc419c063f27");katex.render("T", 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 temperature
$var element = document.getElementById("moose-equation-9369e1e8-55b5-47c2-bf5d-e37fd8d87527");katex.render("T_{\\mathrm{ref}}", 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 reference temperature, which is a coupled variable (usually linear-Lagrange) or a real number
$var element = document.getElementById("moose-equation-bc772a4c-fc29-4047-b0a7-c8bb967de899");katex.render("\\epsilon^{\\mathrm{total}}_{ii}", 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 total volumetric strain
$var element = document.getElementById("moose-equation-04e652b4-8468-4765-9cc0-63306276bc9e");katex.render("\\epsilon^{\\mathrm{total}}_{ii\\ \\mathrm{ref}}", 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 reference total volumetric strain, which is a constant monomial coupled variable or a real number.

In addition, it is possible to place a lower bound on porosity: $var element = document.getElementById("moose-equation-86c9a0be-b119-4954-b364-8e9b59618ea5");katex.render("\\phi \\geq \\phi_{\\mathrm{min}}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . This can be physically useful, but can also be numerically important because during Newton-Raphson iterations, the linear relationship can sometimes yield small, or even negative, values of $var element = document.getElementById("moose-equation-48b15be4-70e6-48e1-81fd-a0ab66b2a136");katex.render("\\phi", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .

Descriptions of other porosity classes can be found in the porosity documentation

Input Parameters

PorousFlowDictatorThe UserObject that holds the list of PorousFlow variable names
C++ Type:UserObjectName
Unit:(no unit assumed)
Controllable:No
Description:The UserObject that holds the list of PorousFlow variable names
porosity_refThe porosity at reference effective porepressure, temperature and volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The porosity at reference effective porepressure, temperature and volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)

Required Parameters

P_coeff0Effective porepressure coefficient
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Effective porepressure coefficient
P_refThe reference effective porepressure. This should usually be a linear-lagrange variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The reference effective porepressure. This should usually be a linear-lagrange variable
T_coeff0Temperature coefficient
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature coefficient
T_refThe reference temperature. This should usually be a linear-lagrange variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The reference temperature. This should usually be a linear-lagrange variable
at_nodesFalseEvaluate Material properties at nodes instead of quadpoints
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Evaluate Material properties at nodes instead of quadpoints
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Unit:(no unit assumed)
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Unit:(no unit assumed)
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object 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
Unit:(no unit assumed)
Controllable:No
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
Unit:(no unit assumed)
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
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
declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
epv_coeff0Volumetric-strain coefficient
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Volumetric-strain coefficient
epv_refThe reference volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The reference volumetric strain. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)
porosity_min0Minimum allowed value of the porosity: if the linear relationship gives values less than this value, then porosity is set to this value instead
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Minimum allowed value of the porosity: if the linear relationship gives values less than this value, then porosity is set to this value instead
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
Unit:(no unit assumed)
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.
strain_at_nearest_qpFalseWhen calculating nodal porosity that depends on strain, use the strain at the nearest quadpoint. This adds a small extra computational burden, and is not necessary for simulations involving only linear lagrange elements. If you set this to true, you will also want to set the same parameter to true for related Kernels and Materials
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:When calculating nodal porosity that depends on strain, use the strain at the nearest quadpoint. This adds a small extra computational burden, and is not necessary for simulations involving only linear lagrange elements. If you set this to true, you will also want to set the same parameter to true for related Kernels and Materials
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
Unit:(no unit assumed)
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.

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Unit:(no unit assumed)
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
Unit:(no unit assumed)
Controllable:No
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
Unit:(no unit assumed)
Controllable:No
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
Unit:(no unit assumed)
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.
zero_modifier0.001If the linear relationship produces porosity < porosity_min, then porosity is set porosity_min. This means the derivatives of it will be zero. However, these zero derivatives often result in poor NR convergence, so the derivatives are set to _zero_modifier * (values that are relevant for min porosity) to hint to the NR that porosity is not always constant
Default:0.001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:If the linear relationship produces porosity < porosity_min, then porosity is set porosity_min. This means the derivatives of it will be zero. However, these zero derivatives often result in poor NR convergence, so the derivatives are set to _zero_modifier * (values that are relevant for min porosity) to hint to the NR that porosity is not always constant

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>
Unit:(no unit assumed)
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Unit:(no unit assumed)
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

Input Files

(modules/porous_flow/test/tests/poroperm/linear_test_vals.i)
(modules/porous_flow/test/tests/poroperm/linear_except1.i)
(modules/porous_flow/examples/multiapp_fracture_flow/3dFracture/fracture_only_aperture_changing.i)
(modules/porous_flow/test/tests/jacobian/linear_por.i)

(modules/porous_flow/test/tests/poroperm/linear_test_vals.i)

# Testing PorousFlowPorosityLinear produces correct values:
# porosity = porosity_ref + P_coeff * (P - P_ref) + T_coeff * (T - T_ref) + epv_coeff * (epv - epv_coeff)
#          = 0.5 + 2 * (1 - 0.5) + 0.5 * (2 - -3) + 4 * (3 - 2.5)
#          = 6
[GlobalParams]
  PorousFlowDictator = dictator
[]

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    number_fluid_phases = 1
    number_fluid_components = 1
    porous_flow_vars = pp
  []
[]

[Variables]
  [pp]
    initial_condition = 1
  []
  [T]
    initial_condition = 2
  []
  [disp]
  []
[]

[ICs]
  [disp]
    type = FunctionIC
    variable = disp
    function = '3 * x'
  []
[]

[Kernels]
  [pp]
    type = TimeDerivative
    variable = pp
  []
  [T]
    type = TimeDerivative
    variable = T
  []
  [disp]
    type = TimeDerivative
    variable = disp
  []
[]

[AuxVariables]
  [porosity]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
  []
[]

[Postprocessors]
  [porosity]
    type = PointValue
    point = '0 0 0'
    variable = porosity
  []
[]

[Materials]
  [ps]
    type = PorousFlow1PhaseFullySaturated
    porepressure = pp
  []
  [temperature]
    type = PorousFlowTemperature
    temperature = T
  []
  [pf]
    type = PorousFlowEffectiveFluidPressure
  []
  [total_strain]
    type = ComputeSmallStrain
    displacements = disp
  []
  [volstrain]
    type = PorousFlowVolumetricStrain
    displacements = disp
  []
  [porosity]
    type = PorousFlowPorosityLinear
    porosity_ref = 0.5
    P_ref = 0.5
    P_coeff = 2.0
    T_ref = -3.0
    T_coeff = 0.5
    epv_ref = 2.5
    epv_coeff = 4.0
  []
[]

[Executioner]
  type = Transient
  dt = 1
  num_steps = 1
[]

[Outputs]
  csv = true
[]

(modules/porous_flow/test/tests/poroperm/linear_except1.i)

# Exception testing of PorousFlowPorosityLinear: demonstrating that an error is thrown if there are missing Materials
[GlobalParams]
  PorousFlowDictator = dictator
[]

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
  []
[]

[Problem]
  kernel_coverage_check = false
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    number_fluid_phases = 1
    number_fluid_components = 1
    porous_flow_vars = pp
  []
[]

[Variables]
  [pp]
  []
  [disp]
  []
[]

[Kernels]
  [pp]
    type = Diffusion
    variable = pp
  []
[]

[Materials]
  [ps]
    type = PorousFlow1PhaseFullySaturated
    porepressure = pp
  []
  [temperature]
    type = PorousFlowTemperature
  []
  [pf]
    type = PorousFlowEffectiveFluidPressure
  []
  [volstrain]
    type = PorousFlowVolumetricStrain
    displacements = pp
  []
  [porosity]
    type = PorousFlowPorosityLinear
    porosity_ref = 0.1
  []
  [total_strain]
    type = ComputeSmallStrain
    displacements = disp
  []
[]

[Executioner]
  type = Transient
  dt = 1
  num_steps = 1
[]

(modules/porous_flow/examples/multiapp_fracture_flow/3dFracture/fracture_only_aperture_changing.i)

# Cold water injection into one side of the fracture network, and production from the other side
injection_rate = 10 # kg/s

[Mesh]
  uniform_refine = 0
  [cluster34]
    type = FileMeshGenerator
    file = 'Cluster_34.exo'
  []
  [injection_node]
    type = BoundingBoxNodeSetGenerator
    input = cluster34
    bottom_left = '-1000 0 -1000'
    top_right = '1000 0.504 1000'
    new_boundary = injection_node
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 -9.81E-6' # Note the value, because of pressure_unit
[]

[Variables]
  [frac_P]
    scaling = 1E6
  []
  [frac_T]
    initial_condition = 473
  []
[]

[ICs]
  [frac_P]
    type = FunctionIC
    variable = frac_P
    function = insitu_pp
  []
[]

[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = frac_P
  temperature = frac_T
  fp = water
  pressure_unit = MPa
[]

[Kernels]
  [toMatrix]
    type = PorousFlowHeatMassTransfer
    variable = frac_T
    v = transferred_matrix_T
    transfer_coefficient = heat_transfer_coefficient
    save_in = joules_per_s
  []
[]

[AuxVariables]
  [heat_transfer_coefficient]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = 0.0
  []
  [transferred_matrix_T]
    initial_condition = 473
  []
  [joules_per_s]
  []
  [normal_dirn_x]
    family = MONOMIAL
    order = CONSTANT
  []
  [normal_dirn_y]
    family = MONOMIAL
    order = CONSTANT
  []
  [normal_dirn_z]
    family = MONOMIAL
    order = CONSTANT
  []
  [enclosing_element_normal_length]
    family = MONOMIAL
    order = CONSTANT
  []
  [enclosing_element_normal_thermal_cond]
    family = MONOMIAL
    order = CONSTANT
  []
  [aperture]
    family = MONOMIAL
    order = CONSTANT
  []
  [perm_times_app]
    family = MONOMIAL
    order = CONSTANT
  []
  [density]
    family = MONOMIAL
    order = CONSTANT
  []
  [viscosity]
    family = MONOMIAL
    order = CONSTANT
  []
  [insitu_pp]
  []
[]

[AuxKernels]
  [normal_dirn_x_auxk]
    type = PorousFlowElementNormal
    variable = normal_dirn_x
    component = x
  []
  [normal_dirn_y]
    type = PorousFlowElementNormal
    variable = normal_dirn_y
    component = y
  []
  [normal_dirn_z]
    type = PorousFlowElementNormal
    variable = normal_dirn_z
    component = z
  []
  [heat_transfer_coefficient_auxk]
    type = ParsedAux
    variable = heat_transfer_coefficient
    coupled_variables = 'enclosing_element_normal_length enclosing_element_normal_thermal_cond'
    constant_names = h_s
    constant_expressions = 1E3 # should be much bigger than thermal_conductivity / L ~ 1
    expression = 'if(enclosing_element_normal_length = 0, 0, h_s * enclosing_element_normal_thermal_cond * 2 * enclosing_element_normal_length / (h_s * enclosing_element_normal_length * enclosing_element_normal_length + enclosing_element_normal_thermal_cond * 2 * enclosing_element_normal_length))'
  []
  [aperture]
    type = PorousFlowPropertyAux
    variable = aperture
    property = porosity
  []
  [perm_times_app]
    type = PorousFlowPropertyAux
    variable = perm_times_app
    property = permeability
    row = 0
    column = 0
  []
  [density]
    type = PorousFlowPropertyAux
    variable = density
    property = density
    phase = 0
  []
  [viscosity]
    type = PorousFlowPropertyAux
    variable = viscosity
    property = viscosity
    phase = 0
  []
  [insitu_pp]
    type = FunctionAux
    execute_on = initial
    variable = insitu_pp
    function = insitu_pp
  []
[]

[BCs]
  [inject_heat]
    type = DirichletBC
    boundary = injection_node
    variable = frac_T
    value = 373
  []
[]

[DiracKernels]
  [inject_fluid]
    type = PorousFlowPointSourceFromPostprocessor
    mass_flux = ${injection_rate}
    point = '58.8124 0.50384 74.7838'
    variable = frac_P
  []
  [withdraw_fluid]
    type = PorousFlowPeacemanBorehole
    SumQuantityUO = kg_out_uo
    bottom_p_or_t = 10.6 # 1MPa + approx insitu at production point, to prevent aperture closing due to low porepressures
    character = 1
    line_length = 1
    point_file = production.xyz
    unit_weight = '0 0 0'
    fluid_phase = 0
    use_mobility = true
    variable = frac_P
  []
  [withdraw_heat]
    type = PorousFlowPeacemanBorehole
    SumQuantityUO = J_out_uo
    bottom_p_or_t = 10.6 # 1MPa + approx insitu at production point, to prevent aperture closing due to low porepressures
    character = 1
    line_length = 1
    point_file = production.xyz
    unit_weight = '0 0 0'
    fluid_phase = 0
    use_mobility = true
    use_enthalpy = true
    variable = frac_T
  []
[]

[UserObjects]
  [kg_out_uo]
    type = PorousFlowSumQuantity
  []
  [J_out_uo]
    type = PorousFlowSumQuantity
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [water]
    type = TabulatedBicubicFluidProperties
    fp = true_water
    temperature_min = 275 # K
    temperature_max = 600
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated.csv
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosityLinear
    porosity_ref = 1E-4 # fracture porosity = 1.0, but must include fracture aperture of 1E-4 at P = insitu_pp
    P_ref = insitu_pp
    P_coeff = 1E-3 # this is in metres/MPa, ie for P_ref = 1/P_coeff, the aperture becomes 1 metre
    porosity_min = 1E-5
  []
  [permeability]
    type = PorousFlowPermeabilityKozenyCarman
    k0 = 1E-15 # fracture perm = 1E-11 m^2, but must include fracture aperture of 1E-4
    poroperm_function = kozeny_carman_phi0
    m = 0
    n = 3
    phi0 = 1E-4
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2700 # kg/m^3
    specific_heat_capacity = 0 # basically no rock inside the fracture
  []
  [aq_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0.6E-4 0 0  0 0.6E-4 0  0 0 0.6E-4' # thermal conductivity of water times fracture aperture.  This should increase linearly with aperture, but is set constant in this model
  []
[]

[Functions]
  [kg_rate]
    type = ParsedFunction
    symbol_values = 'dt kg_out'
    symbol_names = 'dt kg_out'
    expression = 'kg_out/dt'
  []
  [insitu_pp]
    type = ParsedFunction
    expression = '10 - 0.847E-2 * z' # Approximate hydrostatic in MPa
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
    outputs = 'none'
  []
  [kg_out]
    type = PorousFlowPlotQuantity
    uo = kg_out_uo
  []
  [kg_per_s]
    type = FunctionValuePostprocessor
    function = kg_rate
  []
  [J_out]
    type = PorousFlowPlotQuantity
    uo = J_out_uo
  []
  [TK_out]
    type = PointValue
    variable = frac_T
    point = '101.705 160.459 39.5722'
  []
  [P_out]
    type = PointValue
    variable = frac_P
    point = '101.705 160.459 39.5722'
  []
  [P_in]
    type = PointValue
    variable = frac_P
    point = '58.8124 0.50384 74.7838'
  []
[]

[VectorPostprocessors]
  [heat_transfer_rate]
    type = NodalValueSampler
    outputs = none
    sort_by = id
    variable = joules_per_s
  []
[]

[Preconditioning]
  [entire_jacobian]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2             '
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    optimal_iterations = 10
    growth_factor = 1.5
  []
  dtmax = 1E8
  end_time = 1E8
  nl_abs_tol = 1E-3
  nl_max_its = 20
[]

[Outputs]
  print_linear_residuals = false
  csv = true
  [ex]
    type = Exodus
    sync_times = '1 10 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 6100 6200 6300 6400 6500 6600 6700 6800 6900 7000 7100 7200 7300 7400 7500 7600 7700 7800 7900 8000 8100 8200 8300 8400 8500 8600 8700 8800 8900 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 30000 50000 70000 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000 2100000 2200000 2300000 2400000 2500000 2600000 2700000 2800000 2900000'
    sync_only = true
  []
[]

(modules/porous_flow/test/tests/jacobian/linear_por.i)

# Testing Jacobian resulting from PorousFlowPorosityLinear in a THM situation
[GlobalParams]
  PorousFlowDictator = dictator
  strain_at_nearest_qp = true
[]

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Variables]
  [pp]
    initial_condition = 1
  []
  [T]
    initial_condition = 2
  []
  [disp]
  []
[]

[ICs]
  [disp]
    type = FunctionIC
    variable = disp
    function = '3 * x'
  []
[]

[BCs]
  [disp]
    type = FunctionDirichletBC
    boundary = 'left right top bottom front back'
    variable = disp
    function = '3 * x'
  []
[]

[PorousFlowFullySaturated]
  coupling_type = ThermoHydroMechanical
  fp = simple_fluid
  porepressure = pp
  temperature = T
  displacements = 'disp disp disp'
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosityLinear
    porosity_ref = 0.5
    P_ref = 0.5
    P_coeff = 1.0
    T_ref = -3.0
    T_coeff = 1.0
    epv_ref = 2.5
    epv_coeff = 1.0
  []
  [perm]
    type = PorousFlowPermeabilityConst
    permeability = '0 0 0  0 0 0  0 0 0'
  []
  [matrix_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 0.0
    specific_heat_capacity = 0.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0 0 0  0 0 0  0 0 0'
  []
  [density]
    type = GenericConstantMaterial
    prop_names = density
    prop_values = 0.0
  []
  [elasticity]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1E-99
    poissons_ratio = 0
  []
  [strain]
    type = ComputeSmallStrain
    displacements = 'disp disp disp'
  []
  [stress]
    type = ComputeLinearElasticStress
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  num_steps = 1
#  petsc_options = '-snes_test_jacobian -snes_force_iteration'
#  petsc_options_iname = '-snes_type --ksp_type -pc_type -snes_convergence_test'
#  petsc_options_value = ' ksponly     preonly   none     skip'
[]

Input Parameters
Input Files