- allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).Default:False C++ Type:bool Controllable:No Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable). 
- execute_onTIMESTEP_ENDThe list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.Default:TIMESTEP_END C++ Type:ExecFlagEnum Controllable:No Description:The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html. 
- execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.Default:0 C++ Type:int Controllable:No Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group. 
- force_postauxFalseForces the UserObject to be executed in POSTAUXDefault:False C++ Type:bool Controllable:No Description:Forces the UserObject to be executed in POSTAUX 
- force_preauxFalseForces the UserObject to be executed in PREAUXDefault:False C++ Type:bool Controllable:No Description:Forces the UserObject to be executed in PREAUX 
- force_preicFalseForces the UserObject to be executed in PREIC during initial setupDefault:False C++ Type:bool Controllable:No Description:Forces the UserObject to be executed in PREIC during initial setup 
PorousFlowSumQuantity
This is used to record the total fluid (kg) or heat (J) flux that is produced by a PorousFlow DiracKernel in a time step.  See polyline sinks for an extended discussion.
Input 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. 
- enableTrueSet the enabled status of the MooseObject.Default:True C++ Type:bool Controllable:Yes Description:Set the enabled status of the MooseObject. 
- 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 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
- 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. 
- 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. 
Material Property Retrieval Parameters
Input Files
- (modules/porous_flow/test/tests/jacobian/line_sink03.i)
- (modules/combined/examples/geochem-porous_flow/forge/porous_flow.i)
- (modules/porous_flow/test/tests/dirackernels/bh07.i)
- (modules/porous_flow/examples/multiapp_fracture_flow/3dFracture/fracture_only_aperture_changing.i)
- (modules/porous_flow/examples/groundwater/ex01.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except03.i)
- (modules/porous_flow/examples/groundwater/ex02_steady_state.i)
- (modules/porous_flow/test/tests/jacobian/line_sink04.i)
- (modules/porous_flow/test/tests/dirackernels/pls02.i)
- (modules/porous_flow/test/tests/dirackernels/pls02reporter.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except09.i)
- (modules/porous_flow/test/tests/jacobian/line_sink01.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except08.i)
- (modules/porous_flow/test/tests/dirackernels/bh02.i)
- (modules/porous_flow/examples/groundwater/ex02_abstraction.i)
- (modules/porous_flow/test/tests/jacobian/line_sink02.i)
- (modules/porous_flow/test/tests/dirackernels/pls03.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except13.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except04.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except16.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except10.i)
- (modules/porous_flow/test/tests/dirackernels/pls03_action.i)
- (modules/porous_flow/test/tests/dirackernels/strain_at_nearest_qp.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except11.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except02.i)
- (modules/combined/examples/geochem-porous_flow/geotes_2D/porous_flow.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except05.i)
- (modules/combined/examples/geochem-porous_flow/geotes_weber_tensleep/porous_flow.i)
- (modules/porous_flow/test/tests/dirackernels/injection_production.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except06.i)
- (modules/porous_flow/test/tests/dirackernels/pls01.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except01.i)
- (modules/porous_flow/test/tests/dirackernels/bh04.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except14.i)
- (modules/porous_flow/test/tests/actions/basicthm_borehole.i)
- (modules/porous_flow/test/tests/actions/fullsat_borehole.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except12.i)
- (modules/porous_flow/test/tests/dirackernels/bh05.i)
- (modules/porous_flow/test/tests/dirackernels/bh03.i)
- (modules/porous_flow/test/tests/dirackernels/bh02reporter.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except15.i)
- (modules/porous_flow/test/tests/dirackernels/injection_with_plasticity.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except07.i)
(modules/porous_flow/test/tests/jacobian/line_sink03.i)
# PorousFlowPeacemanBorehole with 2-phase, 3-components, with enthalpy, internal_energy, and thermal_conductivity
# NOTE: this test has suffered from repeated failures since its inception.  The problem always appears to be caused by having too many Dirac points in an element: see #10471.  As of Nov2020, the dirac7 DiracKernel uses only one Dirac point, not ten_points.bh.  One day it would be good to be able to use point_file = ten_points.bh
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [ppwater]
  []
  [ppgas]
  []
  [massfrac_ph0_sp0]
  []
  [massfrac_ph0_sp1]
  []
  [massfrac_ph1_sp0]
  []
  [massfrac_ph1_sp1]
  []
  [temp]
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp ppwater ppgas massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
    number_fluid_phases = 2
    number_fluid_components = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
  [dummy_outflow0]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow1]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow2]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow3]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow4]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow5]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow6]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow7]
    type = PorousFlowSumQuantity
  []
[]
[ICs]
  [temp]
    type = RandomIC
    variable = temp
    min = 1
    max = 2
  []
  [ppwater]
    type = RandomIC
    variable = ppwater
    min = -1
    max = 0
  []
  [ppgas]
    type = RandomIC
    variable = ppgas
    min = 0
    max = 1
  []
  [massfrac_ph0_sp0]
    type = RandomIC
    variable = massfrac_ph0_sp0
    min = 0
    max = 1
  []
  [massfrac_ph0_sp1]
    type = RandomIC
    variable = massfrac_ph0_sp1
    min = 0
    max = 1
  []
  [massfrac_ph1_sp0]
    type = RandomIC
    variable = massfrac_ph1_sp0
    min = 0
    max = 1
  []
  [massfrac_ph1_sp1]
    type = RandomIC
    variable = massfrac_ph1_sp1
    min = 0
    max = 1
  []
[]
[Kernels]
  [dummy_temp]
    type = TimeDerivative
    variable = temp
  []
  [dummy_ppwater]
    type = TimeDerivative
    variable = ppwater
  []
  [dummy_ppgas]
    type = TimeDerivative
    variable = ppgas
  []
  [dummy_m00]
    type = TimeDerivative
    variable = massfrac_ph0_sp0
  []
  [dummy_m01]
    type = TimeDerivative
    variable = massfrac_ph0_sp1
  []
  [dummy_m10]
    type = TimeDerivative
    variable = massfrac_ph1_sp0
  []
  [dummy_m11]
    type = TimeDerivative
    variable = massfrac_ph1_sp1
  []
[]
[FluidProperties]
  [simple_fluid0]
    type = SimpleFluidProperties
    bulk_modulus = 1.5
    density0 = 1
    thermal_expansion = 0
    viscosity = 1
    cv = 1.1
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 0.5
    density0 = 0.5
    thermal_expansion = 0
    viscosity = 1.4
    cv = 1.8
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow2PhasePP
    phase0_porepressure = ppwater
    phase1_porepressure = ppgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
  []
  [simple_fluid0]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1 0 0 0 2 0 0 0 3'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0.1 0.02 0.03 0.02 0.0 0.01 0.03 0.01 0.3'
  []
[]
[DiracKernels]
  #active = 'dirac6 dirac2' # incorrect jacobian for ny=2
  #active = 'dirac0 dirac1 dirac2 dirac3 dirac4 dirac5' # correct jacobian for ny=2
  #active = 'dirac0 dirac1 dirac2 dirac3 dirac4 dirac5 dirac6' # incorrect jacobian for ny=2
  #active = 'dirac0 dirac1 dirac2 dirac3 dirac4 dirac5 dirac7' # correct jacobian in dbg, but not in opt, for ny=2
  #active = 'dirac0 dirac1 dirac2 dirac3 dirac4 dirac5 dirac6' # incorrect jacobian in dbg, but correct for opt, for ny=1
  #active = 'dirac0 dirac1 dirac2 dirac3 dirac4 dirac5' # correct jacobian, for ny=1
  #active = 'dirac0 dirac1 dirac2 dirac3 dirac4 dirac5 dirac6' # incorrect jacobian in dbg, but correct for opt, for ny=1.  row24, col 21 and 22 are wrong.  row24=node3, 21=ppwater, 22=ppgas, 24=massfrac_ph0_sp1 (all at node3)
  [dirac0]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = ppwater
    point_file = one_point.bh
    line_length = 1
    SumQuantityUO = dummy_outflow0
    character = 1
    bottom_p_or_t = -10
    unit_weight = '1 2 3'
    re_constant = 0.123
  []
  [dirac1]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = ppgas
    line_length = 1
    line_direction = '-1 -1 -1'
    use_relative_permeability = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow1
    character = -0.5
    bottom_p_or_t = 10
    unit_weight = '1 2 -3'
    re_constant = 0.3
  []
  [dirac2]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = massfrac_ph0_sp0
    line_length = 1.3
    line_direction = '1 0 1'
    use_mobility = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow2
    character = 0.6
    bottom_p_or_t = -4
    unit_weight = '-1 -2 -3'
    re_constant = 0.4
  []
  [dirac3]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = massfrac_ph0_sp1
    line_length = 1.3
    line_direction = '1 1 1'
    use_enthalpy = true
    mass_fraction_component = 0
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow3
    character = -1
    bottom_p_or_t = 3
    unit_weight = '0.1 0.2 0.3'
    re_constant = 0.5
  []
  [dirac4]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = massfrac_ph1_sp0
    function_of = temperature
    line_length = 0.9
    line_direction = '1 1 1'
    mass_fraction_component = 1
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow4
    character = 1.1
    bottom_p_or_t = -7
    unit_weight = '-1 2 3'
    re_constant = 0.6
  []
  [dirac5]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = temp
    line_length = 0.9
    function_of = temperature
    line_direction = '1 2 3'
    mass_fraction_component = 2
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow5
    character = 0.9
    bottom_p_or_t = -8
    unit_weight = '1 2 1'
    re_constant = 0.7
  []
  [dirac6]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = ppwater
    point_file = nine_points.bh
    SumQuantityUO = dummy_outflow6
    character = 0
    bottom_p_or_t = 10
    unit_weight = '0.0 0.0 0.0'
  []
  [dirac7]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = massfrac_ph0_sp0
    use_mobility = true
    mass_fraction_component = 1
    use_relative_permeability = true
    use_internal_energy = true
    point_file = one_point.bh
    #NOTE this commented-out line: point_file = ten_points.bh
    SumQuantityUO = dummy_outflow7
    character = -1
    bottom_p_or_t = 10
    unit_weight = '0.1 0.2 0.3'
  []
[]
[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 = 1
[]
[Outputs]
  file_base = line_sink03
[]
(modules/combined/examples/geochem-porous_flow/forge/porous_flow.i)
# Input file modified from RobPodgorney version
# - 2D instead of 3D with different resolution.  Effectively this means a 1m height of RobPodgorney aquifer is simulated.  RobPodgorney total mass flux is 2.5kg/s meaning 0.25kg/s is appropriate here
# - Celsius instead of Kelvin
# - no use of PorousFlowPointEnthalpySourceFromPostprocessor since that is not yet merged into MOOSE: a DirichletBC is used instead
# - Use of PorousFlowFullySaturated instead of PorousFlowUnsaturated, and the save_component_rate_in feature to record the change in kg of each species at each node for passing to the Geochem simulation
# - MultiApps and Transfers to transfer information between this simulation and the aquifer_geochemistry.i simulation
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 225
    ny = 200
    xmin = -400
    xmax = 500
    ymin = -400
    ymax = 400
  []
  [injection_node]
    input = gen
    type = ExtraNodesetGenerator
    new_boundary = injection_node
    coord = '0 0 0'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]
[Variables]
  [f_H]
    initial_condition = 8.201229858451E-07
  []
  [f_Na]
    initial_condition = 2.281094143525E-03
  []
  [f_K]
    initial_condition = 2.305489507836E-04
  []
  [f_Ca]
    initial_condition = 5.818776782059E-04
  []
  [f_Mg]
    initial_condition = 1.539513498238E-07
  []
  [f_SiO2]
    initial_condition = 2.691822196469E-04
  []
  [f_Al]
    initial_condition = 4.457519474122E-08
  []
  [f_Cl]
    initial_condition = 4.744309776594E-03
  []
  [f_SO4]
    initial_condition = 9.516650880811E-06
  []
  [f_HCO3]
    initial_condition = 5.906126982324E-05
  []
  [porepressure]
    initial_condition = 20E6
  []
  [temperature]
    initial_condition = 220 # degC
    scaling = 1E-6 # fluid enthalpy is roughly 1E6
  []
[]
[BCs]
  [source_temperature]
    type = DirichletBC
    boundary = injection_node
    variable = temperature
    value = 70 # degC
  []
[]
[DiracKernels]
  [inject_H]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 4.790385871045E-08
    variable = f_H
  []
  [inject_Na]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 7.586252963780E-07
    variable = f_Na
  []
  [inject_K]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.746517625125E-07
    variable = f_K
  []
  [inject_Ca]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 7.775129478597E-07
    variable = f_Ca
  []
  [inject_Mg]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 1.749872109005E-07
    variable = f_Mg
  []
  [inject_SiO2]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 4.100547515915E-06
    variable = f_SiO2
  []
  [inject_Al]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.502408592080E-08
    variable = f_Al
  []
  [inject_Cl]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 1.997260386272E-06
    variable = f_Cl
  []
  [inject_SO4]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.497372164191E-07
    variable = f_SO4
  []
  [inject_HCO3]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 5.003150992902E-06
    variable = f_HCO3
  []
  [inject_H2O]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.499865905987E-01
    variable = porepressure
  []
  [produce_H]
    type = PorousFlowPeacemanBorehole
    variable = f_H
    SumQuantityUO = produced_mass_H
    mass_fraction_component = 0
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Na]
    type = PorousFlowPeacemanBorehole
    variable = f_Na
    SumQuantityUO = produced_mass_Na
    mass_fraction_component = 1
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_K]
    type = PorousFlowPeacemanBorehole
    variable = f_K
    SumQuantityUO = produced_mass_K
    mass_fraction_component = 2
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Ca]
    type = PorousFlowPeacemanBorehole
    variable = f_Ca
    SumQuantityUO = produced_mass_Ca
    mass_fraction_component = 3
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Mg]
    type = PorousFlowPeacemanBorehole
    variable = f_Mg
    SumQuantityUO = produced_mass_Mg
    mass_fraction_component = 4
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_SiO2]
    type = PorousFlowPeacemanBorehole
    variable = f_SiO2
    SumQuantityUO = produced_mass_SiO2
    mass_fraction_component = 5
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Al]
    type = PorousFlowPeacemanBorehole
    variable = f_Al
    SumQuantityUO = produced_mass_Al
    mass_fraction_component = 6
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Cl]
    type = PorousFlowPeacemanBorehole
    variable = f_Cl
    SumQuantityUO = produced_mass_Cl
    mass_fraction_component = 7
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_SO4]
    type = PorousFlowPeacemanBorehole
    variable = f_SO4
    SumQuantityUO = produced_mass_SO4
    mass_fraction_component = 8
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_HCO3]
    type = PorousFlowPeacemanBorehole
    variable = f_HCO3
    SumQuantityUO = produced_mass_HCO3
    mass_fraction_component = 9
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_H2O]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = produced_mass_H2O
    mass_fraction_component = 10
    point_file = production.bh
    line_length = 1
    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
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    use_enthalpy = true
    character = 1
  []
[]
[UserObjects]
  [produced_mass_H]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Na]
    type = PorousFlowSumQuantity
  []
  [produced_mass_K]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Ca]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Mg]
    type = PorousFlowSumQuantity
  []
  [produced_mass_SiO2]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Al]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Cl]
    type = PorousFlowSumQuantity
  []
  [produced_mass_SO4]
    type = PorousFlowSumQuantity
  []
  [produced_mass_HCO3]
    type = PorousFlowSumQuantity
  []
  [produced_mass_H2O]
    type = PorousFlowSumQuantity
  []
  [produced_heat]
    type = PorousFlowSumQuantity
  []
[]
[Postprocessors]
  [heat_extracted]
    type = PorousFlowPlotQuantity
    uo = produced_heat
  []
  [approx_production_temperature]
    type = PointValue
    point = '100 0 0'
    variable = temperature
  []
  [mass_extracted_H]
    type = PorousFlowPlotQuantity
    uo = produced_mass_H
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Na]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Na
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_K]
    type = PorousFlowPlotQuantity
    uo = produced_mass_K
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Ca]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Ca
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Mg]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Mg
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_SiO2]
    type = PorousFlowPlotQuantity
    uo = produced_mass_SiO2
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Al]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Al
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Cl]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Cl
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_SO4]
    type = PorousFlowPlotQuantity
    uo = produced_mass_SO4
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_HCO3]
    type = PorousFlowPlotQuantity
    uo = produced_mass_HCO3
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_H2O]
    type = PorousFlowPlotQuantity
    uo = produced_mass_H2O
    execute_on = 'initial timestep_end'
  []
  [mass_extracted]
    type = LinearCombinationPostprocessor
    pp_names = 'mass_extracted_H mass_extracted_Na mass_extracted_K mass_extracted_Ca mass_extracted_Mg mass_extracted_SiO2 mass_extracted_Al mass_extracted_Cl mass_extracted_SO4 mass_extracted_HCO3 mass_extracted_H2O'
    pp_coefs = '1 1 1 1 1 1 1 1 1 1 1'
    execute_on = 'initial timestep_end'
  []
  [dt]
    type = TimestepSize
    execute_on = 'timestep_begin'
  []
[]
[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 2E-4
    bulk_modulus = 2E9
    viscosity = 1E-3
    density0 = 980
    cv = 4000.0
    cp = 4000.0
    porepressure_coefficient = 0
  []
[]
[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = temperature
  mass_fraction_vars = 'f_H f_Na f_K f_Ca f_Mg f_SiO2 f_Al f_Cl f_SO4 f_HCO3'
  save_component_rate_in = 'rate_H rate_Na rate_K rate_Ca rate_Mg rate_SiO2 rate_Al rate_Cl rate_SO4 rate_HCO3 rate_H2O' # change in kg at every node / dt
  fp = the_simple_fluid
  temperature_unit = Celsius
[]
[AuxVariables]
  [rate_H]
  []
  [rate_Na]
  []
  [rate_K]
  []
  [rate_Ca]
  []
  [rate_Mg]
  []
  [rate_SiO2]
  []
  [rate_Al]
  []
  [rate_Cl]
  []
  [rate_SO4]
  []
  [rate_HCO3]
  []
  [rate_H2O]
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.01
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2.5 0 0  0 2.5 0  0 0 2.5'
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    density = 2750.0
    specific_heat_capacity = 900.0
  []
[]
[Preconditioning]
  active = typically_efficient
  [typically_efficient]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_hypre_type'
    petsc_options_value = ' hypre    boomeramg'
  []
  [strong]
    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      ilu           NONZERO                   2'
  []
  [probably_too_strong]
    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 = 31536000 #1 year
  [TimeStepper]
    type = SolutionTimeAdaptiveDT
    dt = 500
  []
[]
[Outputs]
  exodus = true
  csv = true
[]
[MultiApps]
  [react]
    type = TransientMultiApp
    input_files = aquifer_geochemistry.i
    clone_master_mesh = true
    execute_on = 'timestep_end'
  []
[]
[Transfers]
  [changes_due_to_flow]
    type = MultiAppCopyTransfer
    source_variable = 'rate_H rate_Na rate_K rate_Ca rate_Mg rate_SiO2 rate_Al rate_Cl rate_SO4 rate_HCO3 rate_H2O temperature'
    variable = 'pf_rate_H pf_rate_Na pf_rate_K pf_rate_Ca pf_rate_Mg pf_rate_SiO2 pf_rate_Al pf_rate_Cl pf_rate_SO4 pf_rate_HCO3 pf_rate_H2O temperature'
    to_multi_app = react
  []
  [massfrac_from_geochem]
    type = MultiAppCopyTransfer
    source_variable = 'massfrac_H massfrac_Na massfrac_K massfrac_Ca massfrac_Mg massfrac_SiO2 massfrac_Al massfrac_Cl massfrac_SO4 massfrac_HCO3'
    variable = 'f_H f_Na f_K f_Ca f_Mg f_SiO2 f_Al f_Cl f_SO4 f_HCO3'
    from_multi_app = react
  []
[]
(modules/porous_flow/test/tests/dirackernels/bh07.i)
# Comparison with analytical solution for cylindrically-symmetric situation
[Mesh]
  type = FileMesh
  file = bh07_input.e
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Functions]
  [dts]
    type = PiecewiseLinear
    y = '1000 10000'
    x = '100 1000'
  []
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [fflux]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pp
    gravity = '0 0 0'
  []
[]
[BCs]
  [fix_outer]
    type = DirichletBC
    boundary = perimeter
    variable = pp
    value = 1E7
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.8
    alpha = 1e-5
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-11 0 0 0 1E-11 0 0 0 1E-11'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC
    m = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh07.bh
    fluid_phase = 0
    bottom_p_or_t = 0
    unit_weight = '0 0 0'
    use_mobility = true
    re_constant = 0.1594  # use Chen and Zhang version
    character = 2 # double the strength because bh07.bh only fills half the mesh
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
    execute_on = 'initial timestep_end'
  []
  [fluid_mass]
    type = PorousFlowFluidMass
    execute_on = 'initial timestep_end'
  []
[]
[VectorPostprocessors]
  [pp]
    type = LineValueSampler
    variable = pp
    start_point = '0 0 0'
    end_point = '300 0 0'
    sort_by = x
    num_points = 300
    execute_on = timestep_end
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 1E3
  solve_type = NEWTON
  [TimeStepper]
    # get only marginally better results for smaller time steps
    type = FunctionDT
    function = dts
  []
[]
[Outputs]
  file_base = bh07
  [along_line]
    type = CSV
    execute_on = final
  []
  [exodus]
    type = Exodus
    execute_on = 'initial final'
  []
[]
(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/examples/groundwater/ex01.i)
# Groundwater extraction example.
# System consists of two confined aquifers separated by an aquitard
# There is a hydraulic gradient in the upper aquifer
# A well extracts water from the lower aquifer, and the impact on the upper aquifer is observed
# In the center of the model, the roof of the upper aquifer sits 70m below the local water table
[Mesh]
  [basic_mesh]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = -50
    xmax = 50
    nx = 20
    ymin = -25
    ymax = 25
    ny = 10
    zmin = -100
    zmax = -70
    nz = 3
  []
  [lower_aquifer]
    type = SubdomainBoundingBoxGenerator
    input = basic_mesh
    block_id = 1
    block_name = lower_aquifer
    bottom_left = '-1000 -500 -100'
    top_right = '1000 500 -90'
  []
  [aquitard]
    type = SubdomainBoundingBoxGenerator
    input = lower_aquifer
    block_id = 2
    block_name = aquitard
    bottom_left = '-1000 -500 -90'
    top_right = '1000 500 -80'
  []
  [upper_aquifer]
    type = SubdomainBoundingBoxGenerator
    input = aquitard
    block_id = 3
    block_name = upper_aquifer
    bottom_left = '-1000 -500 -80'
    top_right = '1000 500 -70'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
  []
[]
[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = insitu_pp
  []
[]
[BCs]
  [pp]
    type = FunctionDirichletBC
    variable = pp
    function = insitu_pp
    boundary = 'left right top bottom front back'
  []
[]
[Functions]
  [upper_aquifer_head]
    type = ParsedFunction
    expression = '10 + x / 200'
  []
  [lower_aquifer_head]
    type = ParsedFunction
    expression = '20'
  []
  [insitu_head]
    type = ParsedFunction
    symbol_values = 'lower_aquifer_head upper_aquifer_head'
    symbol_names = 'low up'
    expression = 'if(z <= -90, low, if(z >= -80, up, (up * (z + 90) - low * (z + 80)) / (10.0)))'
  []
  [insitu_pp]
    type = ParsedFunction
    symbol_values = 'insitu_head'
    symbol_names = 'h'
    expression = '(h - z) * 1E4'
  []
  [l_rate]
    type = ParsedFunction
    symbol_values = 'm3_produced dt'
    symbol_names = 'm3_produced dt'
    expression = '1000 * m3_produced / dt'
  []
[]
[AuxVariables]
  [insitu_head]
  []
  [head_change]
  []
[]
[AuxKernels]
  [insitu_head]
    type = FunctionAux
    variable = insitu_head
    function = insitu_head
  []
  [head_change]
    type = ParsedAux
    coupled_variables = 'pp insitu_head'
    use_xyzt = true
    expression = 'pp / 1E4 + z - insitu_head'
    variable = head_change
  []
[]
[Postprocessors]
  [m3_produced]
    type = PorousFlowPlotQuantity
    uo = volume_extracted
    outputs = 'none'
  []
  [dt]
    type = TimestepSize
    outputs = 'none'
  []
  [l_per_s]
    type = FunctionValuePostprocessor
    function = l_rate
  []
[]
[VectorPostprocessors]
  [drawdown]
    type = LineValueSampler
    variable = head_change
    start_point = '-50 0 -75'
    end_point = '50 0 -75'
    num_points = 101
    sort_by = x
  []
[]
[PorousFlowBasicTHM]
  fp = simple_fluid
  gravity = '0 0 -10'
  porepressure = pp
  multiply_by_density = false
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    # the following mean that density = 1000 * exp(P / 1E15) ~ 1000
    thermal_expansion = 0
    bulk_modulus = 1E15
  []
[]
[Materials]
  [porosity_aquifers]
    type = PorousFlowPorosityConst
    porosity = 0.05
    block = 'upper_aquifer lower_aquifer'
  []
  [porosity_aquitard]
    type = PorousFlowPorosityConst
    porosity = 0.2
    block = aquitard
  []
  [biot_mod]
    type = PorousFlowConstantBiotModulus
    fluid_bulk_modulus = 2E9
    biot_coefficient = 1.0
  []
  [permeability_aquifers]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
    block = 'upper_aquifer lower_aquifer'
  []
  [permeability_aquitard]
    type = PorousFlowPermeabilityConst
    permeability = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
    block = aquitard
  []
[]
[DiracKernels]
  [sink]
    type = PorousFlowPolyLineSink
    SumQuantityUO = volume_extracted
    point_file = ex01.bh_lower
    line_length = 10
    variable = pp
    # following produces a flux of 0 m^3(water)/m(borehole length)/s if porepressure = 0, and a flux of 1 m^3/m/s if porepressure = 1E9
    p_or_t_vals = '0 1E9'
    fluxes = '0 1'
  []
[]
[UserObjects]
  [volume_extracted]
    type = PorousFlowSumQuantity
  []
[]
[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]
[Executioner]
  type = Transient
  solve_type = Newton
  [TimeStepper]
    type = SolutionTimeAdaptiveDT
    dt = 1.1E5
  []
  end_time = 3.456E5 # 4 days
  nl_abs_tol = 1E-13
[]
[Outputs]
  [csv]
    type = CSV
    file_base = ex01_lower_extraction
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/dirackernels/bh_except03.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
    at_nodes = true # Needed to force expected error
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/examples/groundwater/ex02_steady_state.i)
# Steady-state groundwater model.  See groundwater_models.md for a detailed description
[Mesh]
  [basic_mesh]
    # mesh create by external program: lies within -500<=x<=500 and -200<=y<=200, with varying z
    type = FileMeshGenerator
    file = ex02_mesh.e
  []
  [name_blocks]
    type = RenameBlockGenerator
    input = basic_mesh
    old_block = '2 3 4'
    new_block = 'bot_aquifer aquitard top_aquifer'
  []
  [zmax]
    type = SideSetsFromNormalsGenerator
    input = name_blocks
    normal_tol = 0.1
    new_boundary = zmax
    normals = '0 0 1'
  []
  [xmin_bot_aquifer]
    type = ParsedGenerateSideset
    input = zmax
    included_subdomains = 2
    normal = '-1 0 0'
    combinatorial_geometry = 'x <= -500.0'
    new_sideset_name = xmin_bot_aquifer
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
  []
[]
[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = initial_pp
  []
[]
[BCs]
  [rainfall_recharge]
    type = PorousFlowSink
    boundary = zmax
    variable = pp
    flux_function = -1E-6  # recharge of 0.1mm/day = 1E-4m3/m2/day = 0.1kg/m2/day ~ 1E-6kg/m2/s
  []
  [evapotranspiration]
    type = PorousFlowHalfCubicSink
    boundary = zmax
    variable = pp
    center = 0.0
    cutoff = -5E4 # roots of depth 5m.  5m of water = 5E4 Pa
    use_mobility = true
    fluid_phase = 0
    # Assume pan evaporation of 4mm/day = 4E-3m3/m2/day = 4kg/m2/day ~ 4E-5kg/m2/s
    # Assume that if permeability was 1E-10m^2 and water table at topography then ET acts as pan strength
    # Because use_mobility = true, then 4E-5 = maximum_flux = max * perm * density / visc = max * 1E-4, so max = 40
    max = 40
  []
[]
[DiracKernels]
  [river]
    type = PorousFlowPolyLineSink
    SumQuantityUO = baseflow
    point_file = ex02_river.bh
    # Assume a perennial river.
    # Assume the river has an incision depth of 1m and a stage height of 1.5m, and these are constant in time and uniform over the whole model.  Hence, if groundwater head is 0.5m (5000Pa) there will be no baseflow and leakage.
    p_or_t_vals = '-999995000 5000 1000005000'
    # Assume the riverbed conductance, k_zz*density*river_segment_length*river_width/riverbed_thickness/viscosity = 1E-6*river_segment_length kg/Pa/s
    fluxes = '-1E3 0 1E3'
    variable = pp
  []
[]
[Functions]
  [initial_pp]
    type = SolutionFunction
    scale_factor = 1E4
    from_variable = cosflow_depth
    solution = initial_mesh
  []
  [baseflow_rate]
    type = ParsedFunction
    symbol_names = 'baseflow_kg dt'
    symbol_values = 'baseflow_kg dt'
    expression = 'baseflow_kg / dt * 24.0 * 3600.0 / 400.0'
  []
[]
[PorousFlowUnsaturated]
  fp = simple_fluid
  porepressure = pp
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
  []
[]
[Materials]
  [porosity_everywhere]
    type = PorousFlowPorosityConst
    porosity = 0.05
  []
  [permeability_aquifers]
    type = PorousFlowPermeabilityConst
    block = 'top_aquifer bot_aquifer'
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-13'
  []
  [permeability_aquitard]
    type = PorousFlowPermeabilityConst
    block = aquitard
    permeability = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
  []
[]
[UserObjects]
  [initial_mesh]
    type = SolutionUserObject
    execute_on = INITIAL
    mesh = ex02_mesh.e
    timestep = LATEST
    system_variables = cosflow_depth
  []
  [baseflow]
    type = PorousFlowSumQuantity
  []
[]
[Postprocessors]
  [baseflow_kg]
    type = PorousFlowPlotQuantity
    uo = baseflow
    outputs = 'none'
  []
  [dt]
    type = TimestepSize
    outputs = 'none'
  []
  [baseflow_l_per_m_per_day]
    type = FunctionValuePostprocessor
    function = baseflow_rate
    indirect_dependencies = 'baseflow_kg dt'
  []
[]
[Preconditioning]
  [smp]
    type = SMP
    full = true
    # following 2 lines are not mandatory, but illustrate a popular preconditioner choice in groundwater models
    petsc_options_iname = '-pc_type -sub_pc_type  -pc_asm_overlap'
    petsc_options_value = ' asm      ilu           2              '
  []
[]
[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E6
  [TimeStepper]
    type = FunctionDT
    function = 'max(1E6, t)'
  []
  end_time = 1E12
  nl_abs_tol = 1E-13
[]
[Outputs]
  print_linear_residuals = false
  [ex]
    type = Exodus
    execute_on = final
  []
  [csv]
    type = CSV
  []
[]
(modules/porous_flow/test/tests/jacobian/line_sink04.i)
# PorousFlowPolyLineSink with 2-phase, 3-components, with enthalpy, internal_energy, and thermal_conductivity
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 2
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [ppwater]
  []
  [ppgas]
  []
  [massfrac_ph0_sp0]
  []
  [massfrac_ph0_sp1]
  []
  [massfrac_ph1_sp0]
  []
  [massfrac_ph1_sp1]
  []
  [temp]
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp ppwater ppgas massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
    number_fluid_phases = 2
    number_fluid_components = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
  [dummy_outflow]
    type = PorousFlowSumQuantity
  []
[]
[ICs]
  [temp]
    type = RandomIC
    variable = temp
    min = 1
    max = 2
  []
  [ppwater]
    type = RandomIC
    variable = ppwater
    min = -1
    max = 0
  []
  [ppgas]
    type = RandomIC
    variable = ppgas
    min = 0
    max = 1
  []
  [massfrac_ph0_sp0]
    type = RandomIC
    variable = massfrac_ph0_sp0
    min = 0
    max = 1
  []
  [massfrac_ph0_sp1]
    type = RandomIC
    variable = massfrac_ph0_sp1
    min = 0
    max = 1
  []
  [massfrac_ph1_sp0]
    type = RandomIC
    variable = massfrac_ph1_sp0
    min = 0
    max = 1
  []
  [massfrac_ph1_sp1]
    type = RandomIC
    variable = massfrac_ph1_sp1
    min = 0
    max = 1
  []
[]
[Kernels]
  [dummy_temp]
    type = TimeDerivative
    variable = temp
  []
  [dummy_ppwater]
    type = TimeDerivative
    variable = ppwater
  []
  [dummy_ppgas]
    type = TimeDerivative
    variable = ppgas
  []
  [dummy_m00]
    type = TimeDerivative
    variable = massfrac_ph0_sp0
  []
  [dummy_m01]
    type = TimeDerivative
    variable = massfrac_ph0_sp1
  []
  [dummy_m10]
    type = TimeDerivative
    variable = massfrac_ph1_sp0
  []
  [dummy_m11]
    type = TimeDerivative
    variable = massfrac_ph1_sp1
  []
[]
[FluidProperties]
  [simple_fluid0]
    type = SimpleFluidProperties
    bulk_modulus = 1.5
    density0 = 1
    thermal_expansion = 0
    viscosity = 1
    cv = 1.1
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 0.5
    density0 = 0.5
    thermal_expansion = 0
    viscosity = 1.4
    cv = 1.8
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow2PhasePP
    phase0_porepressure = ppwater
    phase1_porepressure = ppgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
  []
  [simple_fluid0]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1 0 0 0 2 0 0 0 3'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0.1 0.2 0.3 0.2 0 0.1 0.3 0.1 0.1'
  []
[]
[DiracKernels]
  [dirac0]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    variable = ppwater
    point_file = one_point.bh
    line_length = 1
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-0.9 1.5'
    fluxes = '-1.1 2.2'
  []
  [dirac1]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = ppgas
    line_length = 1
    use_relative_permeability = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -2.2'
  []
  [dirac2]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    variable = massfrac_ph0_sp0
    line_length = 1.3
    use_mobility = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac3]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    variable = massfrac_ph0_sp1
    line_length = 1.3
    use_enthalpy = true
    mass_fraction_component = 0
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac4]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = massfrac_ph1_sp0
    function_of = temperature
    line_length = 0.9
    mass_fraction_component = 1
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac5]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = temp
    line_length = 0.9
    mass_fraction_component = 2
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac6]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = massfrac_ph0_sp0
    use_mobility = true
    function_of = temperature
    mass_fraction_component = 1
    use_relative_permeability = true
    use_internal_energy = true
    point_file = ten_points.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '0 -0.2'
  []
[]
[Preconditioning]
  [check]
    type = SMP
    full = true
    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 = 1
[]
[Outputs]
  file_base = line_sink04
[]
(modules/porous_flow/test/tests/dirackernels/pls02.i)
# fully-saturated situation with a poly-line sink with use_mobility=true
# The poly-line consists of 2 points, and has a length
# of 0.5.  Each point is weighted with a weight of 0.1
# The PorousFlowPolyLineSink has
# p_or_t_vals = 0 1E7
# fluxes = 0 1
# so that for 0<=porepressure<=1E7
# base flux = porepressure * 1E-6 * mobility  (measured in kg.m^-1.s^-1),
# and when multiplied by the poly-line length, and
# the weighting of each point, the mass flux is
# flux = porepressure * 0.5*E-8 * mobility (kg.s^-1).
#
# The fluid and matrix properties are:
# porosity = 0.1
# element volume = 8 m^3
# density = dens0 * exp(P / bulk), with bulk = 2E7
# initial porepressure P0 = 1E7
# viscosity = 0.2
# So, fluid mass = 0.8 * density (kg)
#
# The equation to solve is
# d(Mass)/dt = - porepressure * 0.5*E-8 * density / viscosity
#
# PorousFlow discretises time to conserve mass, so to march
# forward in time, we must solve
# Mass(dt) = Mass(0) - P * 0.5E-8 * density / viscosity * dt
# or
# 0.8 * dens0 * exp(P/bulk) = 0.8 * dens0 * exp(P0/bulk) - P * 0.5E-8 * density / viscosity * dt
# For the numbers written above this gives
# P(t=1) = 6.36947 MPa
# which is given precisely by MOOSE
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    viscosity = 0.2
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [pls]
    # This defines a sink that has strength
    # f = L(P) * relperm * L_seg
    # where
    #    L(P) is a piecewise-linear function of porepressure
    #      that is zero at pp=0 and 1 at pp=1E7
    #    relperm is the relative permeability of the fluid
    #    L_seg is the line-segment length associated with
    #      the Dirac points defined in the file pls02.bh
    type = PorousFlowPolyLineSink
    # Because the Variable for this Sink is pp, and pp is associated
    # with the fluid-mass conservation equation, this sink is extracting
    # fluid mass (and not heat energy or something else)
    variable = pp
    # The following specfies that the total fluid mass coming out of
    # the porespace via this sink in this timestep should be recorded
    # in the pls_total_outflow_mass UserObject
    SumQuantityUO = pls_total_outflow_mass
    # The following file defines the polyline geometry
    # which is just two points in this particular example
    point_file = pls02.bh
    # Now define the piecewise-linear function, L
    # First, we want L to be a function of porepressure (and not
    # temperature or something else).  The following means that
    # p_or_t_vals should be intepreted by MOOSE as the zeroth-phase
    # porepressure
    function_of = pressure
    fluid_phase = 0
    # Second, define the piecewise-linear function, L
    # The following means
    #    flux=0 when pp=0  (and also pp<0)
    #    flux=1 when pp=1E7  (and also pp>1E7)
    #    flux=linearly intepolated between pp=0 and pp=1E7
    # When flux>0 this means a sink, while flux<0 means a source
    p_or_t_vals = '0 1E7'
    fluxes = '0 1'
    # Finally, in this case we want to always multiply
    # L by the fluid mobility (of the zeroth phase) and
    # use that in the sink strength instead of the bare L
    # computed above
    use_mobility = true
  []
[]
[Postprocessors]
  [pls_report]
    type = PorousFlowPlotQuantity
    uo = pls_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 pls_report'
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 pls_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]
[Outputs]
  file_base = pls02
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/pls02reporter.i)
# fully-saturated situation with a poly-line sink with use_mobility=true
# The poly-line consists of 2 points, and has a length
# of 0.5.  Each point is weighted with a weight of 0.1
# The PorousFlowPolyLineSink has
# p_or_t_vals = 0 1E7
# fluxes = 0 1
# so that for 0<=porepressure<=1E7
# base flux = porepressure * 1E-6 * mobility  (measured in kg.m^-1.s^-1),
# and when multiplied by the poly-line length, and
# the weighting of each point, the mass flux is
# flux = porepressure * 0.5*E-8 * mobility (kg.s^-1).
#
# The fluid and matrix properties are:
# porosity = 0.1
# element volume = 8 m^3
# density = dens0 * exp(P / bulk), with bulk = 2E7
# initial porepressure P0 = 1E7
# viscosity = 0.2
# So, fluid mass = 0.8 * density (kg)
#
# The equation to solve is
# d(Mass)/dt = - porepressure * 0.5*E-8 * density / viscosity
#
# PorousFlow discretises time to conserve mass, so to march
# forward in time, we must solve
# Mass(dt) = Mass(0) - P * 0.5E-8 * density / viscosity * dt
# or
# 0.8 * dens0 * exp(P/bulk) = 0.8 * dens0 * exp(P0/bulk) - P * 0.5E-8 * density / viscosity * dt
# For the numbers written above this gives
# P(t=1) = 6.36947 MPa
# which is given precisely by MOOSE
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    viscosity = 0.2
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [pls]
    # This defines a sink that has strength
    # f = L(P) * relperm * L_seg
    # where
    #    L(P) is a piecewise-linear function of porepressure
    #      that is zero at pp=0 and 1 at pp=1E7
    #    relperm is the relative permeability of the fluid
    #    L_seg is the line-segment length associated with
    #      the Dirac points defined in the file pls02.bh
    type = PorousFlowPolyLineSink
    # Because the Variable for this Sink is pp, and pp is associated
    # with the fluid-mass conservation equation, this sink is extracting
    # fluid mass (and not heat energy or something else)
    variable = pp
    # The following specfies that the total fluid mass coming out of
    # the porespace via this sink in this timestep should be recorded
    # in the pls_total_outflow_mass UserObject
    SumQuantityUO = pls_total_outflow_mass
    # The following file defines the polyline geometry
    # which is just two points in this particular example
    weight_reporter='pls02file/w'
    x_coord_reporter='pls02file/x'
    y_coord_reporter='pls02file/y'
    z_coord_reporter='pls02file/z'
    # Now define the piecewise-linear function, L
    # First, we want L to be a function of porepressure (and not
    # temperature or something else).  The following means that
    # p_or_t_vals should be intepreted by MOOSE as the zeroth-phase
    # porepressure
    function_of = pressure
    fluid_phase = 0
    # Second, define the piecewise-linear function, L
    # The following means
    #    flux=0 when pp=0  (and also pp<0)
    #    flux=1 when pp=1E7  (and also pp>1E7)
    #    flux=linearly intepolated between pp=0 and pp=1E7
    # When flux>0 this means a sink, while flux<0 means a source
    p_or_t_vals = '0 1E7'
    fluxes = '0 1'
    # Finally, in this case we want to always multiply
    # L by the fluid mobility (of the zeroth phase) and
    # use that in the sink strength instead of the bare L
    # computed above
    use_mobility = true
  []
[]
[Reporters]
  [pls02file]
    # contains contents from pls02.bh
    type=ConstantReporter
    real_vector_names = 'w x y z'
    real_vector_values = '0.10 0.10;
                          0.00 0.00;
                          0.00 0.00;
                         -0.25 0.25'
  []
[]
[Postprocessors]
  [pls_report]
    type = PorousFlowPlotQuantity
    uo = pls_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 pls_report'
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 pls_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]
[Outputs]
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except09.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
    compute_enthalpy = false
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    use_mobility = true
    use_enthalpy = true
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/jacobian/line_sink01.i)
# PorousFlowPeacemanBorehole with 2-phase, 3-components, with enthalpy, internal_energy, and thermal_conductivity
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [ppwater]
  []
  [ppgas]
  []
  [massfrac_ph0_sp0]
  []
  [massfrac_ph0_sp1]
  []
  [massfrac_ph1_sp0]
  []
  [massfrac_ph1_sp1]
  []
  [temp]
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp ppwater ppgas massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
    number_fluid_phases = 2
    number_fluid_components = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
  [dummy_outflow0]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow1]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow2]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow3]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow4]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow5]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow6]
    type = PorousFlowSumQuantity
  []
  [dummy_outflow7]
    type = PorousFlowSumQuantity
  []
[]
[ICs]
  [temp]
    type = RandomIC
    variable = temp
    min = 1
    max = 2
  []
  [ppwater]
    type = RandomIC
    variable = ppwater
    min = -1
    max = 0
  []
  [ppgas]
    type = RandomIC
    variable = ppgas
    min = 0
    max = 1
  []
  [massfrac_ph0_sp0]
    type = RandomIC
    variable = massfrac_ph0_sp0
    min = 0
    max = 1
  []
  [massfrac_ph0_sp1]
    type = RandomIC
    variable = massfrac_ph0_sp1
    min = 0
    max = 1
  []
  [massfrac_ph1_sp0]
    type = RandomIC
    variable = massfrac_ph1_sp0
    min = 0
    max = 1
  []
  [massfrac_ph1_sp1]
    type = RandomIC
    variable = massfrac_ph1_sp1
    min = 0
    max = 1
  []
[]
[Kernels]
  [dummy_temp]
    type = TimeDerivative
    variable = temp
  []
  [dummy_ppwater]
    type = TimeDerivative
    variable = ppwater
  []
  [dummy_ppgas]
    type = TimeDerivative
    variable = ppgas
  []
  [dummy_m00]
    type = TimeDerivative
    variable = massfrac_ph0_sp0
  []
  [dummy_m01]
    type = TimeDerivative
    variable = massfrac_ph0_sp1
  []
  [dummy_m10]
    type = TimeDerivative
    variable = massfrac_ph1_sp0
  []
  [dummy_m11]
    type = TimeDerivative
    variable = massfrac_ph1_sp1
  []
[]
[FluidProperties]
  [simple_fluid0]
    type = SimpleFluidProperties
    bulk_modulus = 1.5
    density0 = 1
    thermal_expansion = 0
    viscosity = 1
    cv = 1.1
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 0.5
    density0 = 0.5
    thermal_expansion = 0
    viscosity = 1.4
    cv = 1.8
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow2PhasePP
    phase0_porepressure = ppwater
    phase1_porepressure = ppgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
  []
  [simple_fluid0]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1 0 0 0 2 0 0 0 3'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0.1 0.02 0.03 0.02 0.0 0.01 0.03 0.01 0.3'
  []
[]
[DiracKernels]
  [dirac0]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = ppwater
    point_file = one_point.bh
    line_length = 1
    SumQuantityUO = dummy_outflow0
    character = 1
    bottom_p_or_t = -10
    unit_weight = '1 2 3'
    re_constant = 0.123
  []
  [dirac1]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = ppgas
    line_length = 1
    line_direction = '-1 -1 -1'
    use_relative_permeability = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow1
    character = -0.5
    bottom_p_or_t = 10
    unit_weight = '1 2 -3'
    re_constant = 0.3
  []
  [dirac2]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = massfrac_ph0_sp0
    line_length = 1.3
    line_direction = '1 0 1'
    use_mobility = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow2
    character = 0.6
    bottom_p_or_t = -4
    unit_weight = '-1 -2 -3'
    re_constant = 0.4
  []
  [dirac3]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = massfrac_ph0_sp1
    line_length = 1.3
    line_direction = '1 1 1'
    use_enthalpy = true
    mass_fraction_component = 0
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow3
    character = -1
    bottom_p_or_t = 3
    unit_weight = '0.1 0.2 0.3'
    re_constant = 0.5
  []
  [dirac4]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = massfrac_ph1_sp0
    function_of = temperature
    line_length = 0.9
    line_direction = '1 1 1'
    mass_fraction_component = 1
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow4
    character = 1.1
    bottom_p_or_t = -7
    unit_weight = '-1 2 3'
    re_constant = 0.6
  []
  [dirac5]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = temp
    line_length = 0.9
    function_of = temperature
    line_direction = '1 2 3'
    mass_fraction_component = 2
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow5
    character = 0.9
    bottom_p_or_t = -8
    unit_weight = '1 2 1'
    re_constant = 0.7
  []
  [dirac6]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 0
    variable = ppwater
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow6
    character = 0
    bottom_p_or_t = 10
    unit_weight = '0.0 0.0 0.0'
  []
  [dirac7]
    type = PorousFlowPeacemanBorehole
    fluid_phase = 1
    variable = massfrac_ph0_sp0
    use_mobility = true
    mass_fraction_component = 1
    use_relative_permeability = true
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow7
    character = -1
    bottom_p_or_t = 10
    unit_weight = '0.1 0.2 0.3'
  []
[]
[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 = 1
[]
[Outputs]
  file_base = line_sink01
[]
(modules/porous_flow/test/tests/dirackernels/bh_except08.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
    at_nodes = false # Needed to force expected error
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    use_mobility = true
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh02.i)
# fully-saturated
# production
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    # Because the Variable for this Sink is pp, and pp is associated
    # with the fluid-mass conservation equation, this sink is extracting
    # fluid mass (and not heat energy or something else)
    variable = pp
    # The following specfies that the total fluid mass coming out of
    # the porespace via this sink in this timestep should be recorded
    # in the pls_total_outflow_mass UserObject
    SumQuantityUO = borehole_total_outflow_mass
    # The following file defines the polyline geometry
    # which is just two points in this particular example
    point_file = bh02.bh
    # First, we want Peacemans f to be a function of porepressure (and not
    # temperature or something else).  So bottom_p_or_t is actually porepressure
    function_of = pressure
    fluid_phase = 0
    # The bottomhole pressure
    bottom_p_or_t = 0
    # In this example there is no increase of the wellbore pressure
    # due to gravity:
    unit_weight = '0 0 0'
    # PeacemanBoreholes should almost always have use_mobility = true
    use_mobility = true
    # This is a production wellbore (a sink of fluid that removes fluid from porespace)
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 bh_report'
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
[Outputs]
  file_base = bh02
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/examples/groundwater/ex02_abstraction.i)
# Abstraction groundwater model.  See groundwater_models.md for a detailed description
[Mesh]
  [from_steady_state]
    type = FileMeshGenerator
    file = gold/ex02_steady_state_ex.e
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
  []
[]
[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = steady_state_pp
  []
[]
[BCs]
  [rainfall_recharge]
    type = PorousFlowSink
    boundary = zmax
    variable = pp
    flux_function = -1E-6  # recharge of 0.1mm/day = 1E-4m3/m2/day = 0.1kg/m2/day ~ 1E-6kg/m2/s
  []
  [evapotranspiration]
    type = PorousFlowHalfCubicSink
    boundary = zmax
    variable = pp
    center = 0.0
    cutoff = -5E4 # roots of depth 5m.  5m of water = 5E4 Pa
    use_mobility = true
    fluid_phase = 0
    # Assume pan evaporation of 4mm/day = 4E-3m3/m2/day = 4kg/m2/day ~ 4E-5kg/m2/s
    # Assume that if permeability was 1E-10m^2 and water table at topography then ET acts as pan strength
    # Because use_mobility = true, then 4E-5 = maximum_flux = max * perm * density / visc = max * 1E-4, so max = 40
    max = 40
  []
[]
[DiracKernels]
  inactive = polyline_sink_borehole
  [river]
    type = PorousFlowPolyLineSink
    SumQuantityUO = baseflow
    point_file = ex02_river.bh
    # Assume a perennial river.
    # Assume the river has an incision depth of 1m and a stage height of 1.5m, and these are constant in time and uniform over the whole model.  Hence, if groundwater head is 0.5m (5000Pa) there will be no baseflow and leakage.
    p_or_t_vals = '-999995000 5000 1000005000'
    # Assume the riverbed conductance, k_zz*density*river_segment_length*river_width/riverbed_thickness/viscosity = 1E-6*river_segment_length kg/Pa/s
    fluxes = '-1E3 0 1E3'
    variable = pp
  []
  [horizontal_borehole]
    type = PorousFlowPeacemanBorehole
    SumQuantityUO = abstraction
    bottom_p_or_t = -1E5
    unit_weight = '0 0 -1E4'
    character = 1.0
    point_file = ex02.bh
    variable = pp
  []
  [polyline_sink_borehole]
    type = PorousFlowPolyLineSink
    SumQuantityUO = abstraction
    fluxes = '-0.4 0 0.4'
    p_or_t_vals = '-1E8 0 1E8'
    point_file = ex02.bh
    variable = pp
  []
[]
[Functions]
  [steady_state_pp]
    type = SolutionFunction
    from_variable = pp
    solution = steady_state_solution
  []
  [baseflow_rate]
    type = ParsedFunction
    symbol_names = 'baseflow_kg dt'
    symbol_values = 'baseflow_kg dt'
    expression = 'baseflow_kg / dt * 24.0 * 3600.0 / 400.0'
  []
  [abstraction_rate]
    type = ParsedFunction
    symbol_names = 'abstraction_kg dt'
    symbol_values = 'abstraction_kg dt'
    expression = 'abstraction_kg / dt * 24.0 * 3600.0'
  []
[]
[AuxVariables]
  [ini_pp]
  []
  [pp_change]
  []
[]
[AuxKernels]
  [ini_pp]
    type = FunctionAux
    variable = ini_pp
    function = steady_state_pp
    execute_on = INITIAL
  []
  [pp_change]
    type = ParsedAux
    variable = pp_change
    coupled_variables = 'pp ini_pp'
    expression = 'pp - ini_pp'
  []
[]
[PorousFlowUnsaturated]
  fp = simple_fluid
  porepressure = pp
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
  []
[]
[Materials]
  [porosity_everywhere]
    type = PorousFlowPorosityConst
    porosity = 0.05
  []
  [permeability_aquifers]
    type = PorousFlowPermeabilityConst
    block = 'top_aquifer bot_aquifer'
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-13'
  []
  [permeability_aquitard]
    type = PorousFlowPermeabilityConst
    block = aquitard
    permeability = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
  []
[]
[UserObjects]
  [steady_state_solution]
    type = SolutionUserObject
    execute_on = INITIAL
    mesh = gold/ex02_steady_state_ex.e
    timestep = LATEST
    system_variables = pp
  []
  [baseflow]
    type = PorousFlowSumQuantity
  []
  [abstraction]
    type = PorousFlowSumQuantity
  []
[]
[Postprocessors]
  [baseflow_kg]
    type = PorousFlowPlotQuantity
    uo = baseflow
    outputs = 'none'
  []
  [dt]
    type = TimestepSize
    outputs = 'none'
  []
  [baseflow_l_per_m_per_day]
    type = FunctionValuePostprocessor
    function = baseflow_rate
    indirect_dependencies = 'baseflow_kg dt'
  []
  [abstraction_kg]
    type = PorousFlowPlotQuantity
    uo = abstraction
    outputs = 'none'
  []
  [abstraction_kg_per_day]
    type = FunctionValuePostprocessor
    function = abstraction_rate
    indirect_dependencies = 'abstraction_kg dt'
  []
[]
[Preconditioning]
  [smp]
    type = SMP
    full = true
    # following 2 lines are not mandatory, but illustrate a popular preconditioner choice in groundwater models
    petsc_options_iname = '-pc_type -sub_pc_type  -pc_asm_overlap'
    petsc_options_value = ' asm      ilu           2              '
  []
[]
[Executioner]
  type = Transient
  solve_type = Newton
  dt = 100
  [TimeStepper]
    type = FunctionDT
    function = 'max(100, t)'
  []
  end_time = 8.64E5 # 10 days
  nl_abs_tol = 1E-11
[]
[Outputs]
  print_linear_residuals = false
  [ex]
    type = Exodus
    execute_on = final
  []
  [csv]
    type = CSV
  []
[]
(modules/porous_flow/test/tests/jacobian/line_sink02.i)
# PorousFlowPolyLineSink with 2-phase, 3-components, with enthalpy, internal_energy, and thermal_conductivity
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 2
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [ppwater]
  []
  [ppgas]
  []
  [massfrac_ph0_sp0]
  []
  [massfrac_ph0_sp1]
  []
  [massfrac_ph1_sp0]
  []
  [massfrac_ph1_sp1]
  []
  [temp]
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp ppwater ppgas massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
    number_fluid_phases = 2
    number_fluid_components = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
  [dummy_outflow]
    type = PorousFlowSumQuantity
  []
[]
[ICs]
  [temp]
    type = RandomIC
    variable = temp
    min = 1
    max = 2
  []
  [ppwater]
    type = RandomIC
    variable = ppwater
    min = -1
    max = 0
  []
  [ppgas]
    type = RandomIC
    variable = ppgas
    min = 0
    max = 1
  []
  [massfrac_ph0_sp0]
    type = RandomIC
    variable = massfrac_ph0_sp0
    min = 0
    max = 1
  []
  [massfrac_ph0_sp1]
    type = RandomIC
    variable = massfrac_ph0_sp1
    min = 0
    max = 1
  []
  [massfrac_ph1_sp0]
    type = RandomIC
    variable = massfrac_ph1_sp0
    min = 0
    max = 1
  []
  [massfrac_ph1_sp1]
    type = RandomIC
    variable = massfrac_ph1_sp1
    min = 0
    max = 1
  []
[]
[Kernels]
  [dummy_temp]
    type = TimeDerivative
    variable = temp
  []
  [dummy_ppwater]
    type = TimeDerivative
    variable = ppwater
  []
  [dummy_ppgas]
    type = TimeDerivative
    variable = ppgas
  []
  [dummy_m00]
    type = TimeDerivative
    variable = massfrac_ph0_sp0
  []
  [dummy_m01]
    type = TimeDerivative
    variable = massfrac_ph0_sp1
  []
  [dummy_m10]
    type = TimeDerivative
    variable = massfrac_ph1_sp0
  []
  [dummy_m11]
    type = TimeDerivative
    variable = massfrac_ph1_sp1
  []
[]
[FluidProperties]
  [simple_fluid0]
    type = SimpleFluidProperties
    bulk_modulus = 1.5
    density0 = 1
    thermal_expansion = 0
    viscosity = 1
    cv = 1.1
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 0.5
    density0 = 0.5
    thermal_expansion = 0
    viscosity = 1.4
    cv = 1.8
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow2PhasePP
    phase0_porepressure = ppwater
    phase1_porepressure = ppgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
  []
  [simple_fluid0]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1 0 0 0 2 0 0 0 3'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0.1 0.2 0.3 0.2 0 0.1 0.3 0.1 0.1'
  []
[]
[DiracKernels]
  [dirac0]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    variable = ppwater
    point_file = one_point.bh
    line_length = 1
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-0.9 1.5'
    fluxes = '-1.1 2.2'
  []
  [dirac1]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = ppgas
    line_length = 1
    use_relative_permeability = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -2.2'
  []
  [dirac2]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    variable = massfrac_ph0_sp0
    line_length = 1.3
    use_mobility = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac3]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    variable = massfrac_ph0_sp1
    line_length = 1.3
    use_enthalpy = true
    mass_fraction_component = 0
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac4]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = massfrac_ph1_sp0
    function_of = temperature
    line_length = 0.9
    mass_fraction_component = 1
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac5]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = temp
    line_length = 0.9
    mass_fraction_component = 2
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '1.1 -0.2'
  []
  [dirac6]
    type = PorousFlowPolyLineSink
    fluid_phase = 1
    variable = massfrac_ph0_sp0
    use_mobility = true
    function_of = temperature
    mass_fraction_component = 1
    use_relative_permeability = true
    use_internal_energy = true
    point_file = one_point.bh
    SumQuantityUO = dummy_outflow
    p_or_t_vals = '-1.9 1.5'
    fluxes = '0 -0.2'
  []
[]
[Preconditioning]
  [check]
    type = SMP
    full = true
    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 = 1
[]
[Outputs]
  file_base = line_sink02
[]
(modules/porous_flow/test/tests/dirackernels/pls03.i)
# Test that the upwinding works correctly.
#
# A poly-line sink sits at the centre of the element.
# It has length=4 and weight=0.5, and extracts fluid
# at a constant rate of
# (1 * relative_permeability) kg.m^-1.s^-1
# Since it sits at the centre of the element, it extracts
# equally from each node, so the rate of extraction from
# each node is
# (0.5 * relative_permeability) kg.s^-1
# including the length and weight effects.
#
# There is no fluid flow.
#
# The initial conditions are such that all nodes have
# relative_permeability=0, except for one which has
# relative_permeaility = 1.  Therefore, all nodes should
# remain at their initial porepressure, except the one.
#
# The porosity is 0.1, and the elemental volume is 2,
# so the fluid mass at the node in question = 0.2 * density / 4,
# where the 4 is the number of nodes in the element.
# In this simulation density = dens0 * exp(P / bulk), with
# dens0 = 100, and bulk = 20 MPa.
# The initial porepressure P0 = 10 MPa, so the final (after
# 1 second of simulation) is
# P(t=1) = 8.748592 MPa
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 0
  xmax = 2
  ymin = 0
  ymax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
  []
[]
[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    #function = if((x<1)&(y<0.5),1E7,-1E7)
    function = if((x<1)&(y>0.5),1E7,-1E7)
    #function = if((x>1)&(y<0.5),1E7,-1E7)
    #function = if((x>1)&(y>0.5),1E7,-1E7)
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    density0 = 100
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC
    phase = 0
    m = 2
    s_res = 0.99
    sum_s_res = 0.99
  []
[]
[DiracKernels]
  [pls]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    point_file = pls03.bh
    use_relative_permeability = true
    line_length = 4
    SumQuantityUO = pls_total_outflow_mass
    variable = pp
    p_or_t_vals = '0 1E7'
    fluxes = '1 1'
  []
[]
[Postprocessors]
  [pls_report]
    type = PorousFlowPlotQuantity
    uo = pls_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 pls_report'
  []
  [p00]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
  [p01]
    type = PointValue
    variable = pp
    point = '0 1 0'
    execute_on = timestep_end
  []
  [p20]
    type = PointValue
    variable = pp
    point = '2 0 0'
    execute_on = timestep_end
  []
  [p21]
    type = PointValue
    variable = pp
    point = '2 1 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 pls_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]
[Outputs]
  file_base = pls03
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except13.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = coincident_points.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except04.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
    at_nodes = true # Needed to force exepected error
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    function_of = temperature
    point_file = bh02.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except16.i)
# fully-saturated
# production
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    function_of = temperature
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except10.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
    compute_internal_energy = false
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    use_mobility = true
    use_internal_energy = true
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/pls03_action.i)
# Test that the upwinding works correctly.
#
# A poly-line sink sits at the centre of the element.
# It has length=4 and weight=0.5, and extracts fluid
# at a constant rate of
# (1 * relative_permeability) kg.m^-1.s^-1
# Since it sits at the centre of the element, it extracts
# equally from each node, so the rate of extraction from
# each node is
# (0.5 * relative_permeability) kg.s^-1
# including the length and weight effects.
#
# There is no fluid flow.
#
# The initial conditions are such that all nodes have
# relative_permeability=0, except for one which has
# relative_permeaility = 1.  Therefore, all nodes should
# remain at their initial porepressure, except the one.
#
# The porosity is 0.1, and the elemental volume is 2,
# so the fluid mass at the node in question = 0.2 * density / 4,
# where the 4 is the number of nodes in the element.
# In this simulation density = dens0 * exp(P / bulk), with
# dens0 = 100, and bulk = 20 MPa.
# The initial porepressure P0 = 10 MPa, so the final (after
# 1 second of simulation) is
# P(t=1) = 8.748592 MPa
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 0
  xmax = 2
  ymin = 0
  ymax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
  []
[]
[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 0.0
    bulk_modulus = 2.0E7
    viscosity = 1.0
    density0 = 100.0
  []
[]
[PorousFlowUnsaturated]
  porepressure = pp
  gravity = '0 0 0'
  fp = the_simple_fluid
  van_genuchten_alpha = 1.0E-7
  van_genuchten_m = 0.5
  relative_permeability_exponent = 2
  residual_saturation = 0.99
[]
[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    #function = if((x<1)&(y<0.5),1E7,-1E7)
    function = if((x<1)&(y>0.5),1E7,-1E7)
    #function = if((x>1)&(y<0.5),1E7,-1E7)
    #function = if((x>1)&(y>0.5),1E7,-1E7)
  []
[]
[UserObjects]
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
[]
[Materials]
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0 0 0  0 0 0  0 0 0'
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]
[DiracKernels]
  [pls]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    point_file = pls03.bh
    use_relative_permeability = true
    line_length = 4
    SumQuantityUO = pls_total_outflow_mass
    variable = pp
    p_or_t_vals = '0 1E7'
    fluxes = '1 1'
  []
[]
[Postprocessors]
  [pls_report]
    type = PorousFlowPlotQuantity
    uo = pls_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 pls_report'
  []
  [p00]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
  [p01]
    type = PointValue
    variable = pp
    point = '0 1 0'
    execute_on = timestep_end
  []
  [p20]
    type = PointValue
    variable = pp
    point = '2 0 0'
    execute_on = timestep_end
  []
  [p21]
    type = PointValue
    variable = pp
    point = '2 1 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 pls_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]
[Outputs]
  file_base = pls03_action
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/strain_at_nearest_qp.i)
# Demonstrates the correct usage of strain_at_nearest_qp when using a nodal PorousFlowPorosity
# For the PorousFlowPorosity Material to require the strain_at_nearest_qp=true flag, it must:
#  - be a nodal Material
#  - be coupled to solid mechanics (mechanical=true)
#  - be part of a simulation with DiracKernels
# The reason for this requirement is that the volumetric strain is a standard Material (at_nodes=false)
# so that it is evaluated at the single Dirac quadpoint, and has size = 1 (assuming just one Dirac point).
# However, the PorousFlowPorosity Material will have size = 2  (number of nodes in the element containing the Dirac point).
# So when the PorousFlowPorosity Material is evaluated, it will use _vol_strain at 2 points.
# If strain_at_nearest_qp=false, then _vol_strain will be evaluated at two quadpoints, but it only has size=1, leading to a segfault
# If strain_at_nearest_qp=true, then _vol_strain will be evaluated correctly just at the single quadpoint
#
# This input file solves no useful physics: it is just illustrating the above point
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
  strain_at_nearest_qp = true
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'disp_x'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [dummy_sum]
    type = PorousFlowSumQuantity
  []
[]
[Variables]
  [disp_x]
  []
[]
[Kernels]
  [dummy]
    type = Diffusion
    variable = disp_x
  []
[]
[DiracKernels]
  [line_sink]
    type = PorousFlowPolyLineSink
    function_of = temperature
    SumQuantityUO = dummy_sum
    point_file = strain_at_nearest_qp.bh
    p_or_t_vals = '0'
    fluxes = '0'
    variable = disp_x
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature  # needed because of the PorousFlowPolyLineSink
  []
  [total_strain]
    type = ComputeSmallStrain
    displacements = disp_x
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
    displacements = disp_x
  []
  [porosity_at_nodes]
    type = PorousFlowPorosity
    mechanical = true # to ensure coupling with volumetric strain
    at_nodes = true  # to ensure evaluation at nodes
    porosity_zero = 0
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
  []
[]
[Executioner]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except11.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    use_relative_permeability = true
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except02.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    mass_fraction_component = 1
    point_file = bh02.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/combined/examples/geochem-porous_flow/geotes_2D/porous_flow.i)
# PorousFlow simulation of injection and production in a 2D aquifer
# Much of this file is standard porous-flow stuff.  The unusual aspects are:
# - transfer of the rates of changes of each species (kg/s) to the aquifer_geochemistry.i simulation.  This is achieved by saving these changes from the PorousFlowMassTimeDerivative residuals
# - transfer of the temperature field to the aquifer_geochemistry.i simulation
# Interesting behaviour can be simulated by this file without its "parent" simulation, exchanger.i.  exchanger.i provides mass-fractions injected via the injection_rate_massfrac_* variables, but since these are more-or-less constant throughout the duration of the exchanger.i simulation, the initial_conditions specified below may be used.  Similar, exchanger.i provides injection_temperature, but that is also constant.
injection_rate = -1.0 # kg/s/m, negative because injection as a source
production_rate = 1.0 # kg/s/m
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 14 # for better resolution, use 56 or 112
    ny = 8  # for better resolution, use 32 or 64
    xmin = -70
    xmax = 70
    ymin = -40
    ymax = 40
  []
  [injection_node]
    input = gen
    type = ExtraNodesetGenerator
    new_boundary = injection_node
    coord = '-30 0 0'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]
[Variables]
  [f0]
    initial_condition = 0.002285946
  []
  [f1]
    initial_condition = 0.0035252
  []
  [f2]
    initial_condition = 1.3741E-05
  []
  [porepressure]
    initial_condition = 2E6
  []
  [temperature]
    initial_condition = 50
    scaling = 1E-6 # fluid enthalpy is roughly 1E6
  []
[]
[BCs]
  [injection_temperature]
    type = MatchedValueBC
    variable = temperature
    v = injection_temperature
    boundary = injection_node
  []
[]
[DiracKernels]
  [inject_Na]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    multiplying_var = injection_rate_massfrac_Na
    point_file = injection.bh
    variable = f0
  []
  [inject_Cl]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    multiplying_var = injection_rate_massfrac_Cl
    point_file = injection.bh
    variable = f1
  []
  [inject_SiO2]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    multiplying_var = injection_rate_massfrac_SiO2
    point_file = injection.bh
    variable = f2
  []
  [inject_H2O]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    multiplying_var = injection_rate_massfrac_H2O
    point_file = injection.bh
    variable = porepressure
  []
  [produce_Na]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Na
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    mass_fraction_component = 0
    point_file = production.bh
    variable = f0
  []
  [produce_Cl]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Cl
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    mass_fraction_component = 1
    point_file = production.bh
    variable = f1
  []
  [produce_SiO2]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_SiO2
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    mass_fraction_component = 2
    point_file = production.bh
    variable = f2
  []
  [produce_H2O]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_H2O
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    mass_fraction_component = 3
    point_file = production.bh
    variable = porepressure
  []
  [produce_heat]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_heat
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    line_length = 1.0
    use_enthalpy = true
    point_file = production.bh
    variable = temperature
  []
[]
[UserObjects]
  [injected_mass]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Na]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Cl]
    type = PorousFlowSumQuantity
  []
  [produced_mass_SiO2]
    type = PorousFlowSumQuantity
  []
  [produced_mass_H2O]
    type = PorousFlowSumQuantity
  []
  [produced_heat]
    type = PorousFlowSumQuantity
  []
[]
[Postprocessors]
  [dt]
    type = TimestepSize
    execute_on = TIMESTEP_BEGIN
  []
  [tot_kg_injected_this_timestep]
    type = PorousFlowPlotQuantity
    uo = injected_mass
  []
  [kg_Na_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Na
  []
  [kg_Cl_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Cl
  []
  [kg_SiO2_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_SiO2
  []
  [kg_H2O_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_H2O
  []
  [mole_rate_Na_produced]
    type = FunctionValuePostprocessor
    function = moles_Na
    indirect_dependencies = 'kg_Na_produced_this_timestep dt'
  []
  [mole_rate_Cl_produced]
    type = FunctionValuePostprocessor
    function = moles_Cl
    indirect_dependencies = 'kg_Cl_produced_this_timestep dt'
  []
  [mole_rate_SiO2_produced]
    type = FunctionValuePostprocessor
    function = moles_SiO2
    indirect_dependencies = 'kg_SiO2_produced_this_timestep dt'
  []
  [mole_rate_H2O_produced]
    type = FunctionValuePostprocessor
    function = moles_H2O
    indirect_dependencies = 'kg_H2O_produced_this_timestep dt'
  []
  [heat_joules_extracted_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_heat
  []
  [production_temperature]
    type = PointValue
    point = '30 0 0'
    variable = temperature
  []
[]
[Functions]
  [moles_Na]
    type = ParsedFunction
    symbol_names = 'kg_Na dt'
    symbol_values = 'kg_Na_produced_this_timestep dt'
    expression = 'kg_Na * 1000 / 22.9898 / dt'
  []
  [moles_Cl]
    type = ParsedFunction
    symbol_names = 'kg_Cl dt'
    symbol_values = 'kg_Cl_produced_this_timestep dt'
    expression = 'kg_Cl * 1000 / 35.453 / dt'
  []
  [moles_SiO2]
    type = ParsedFunction
    symbol_names = 'kg_SiO2 dt'
    symbol_values = 'kg_SiO2_produced_this_timestep dt'
    expression = 'kg_SiO2 * 1000 / 60.0843 / dt'
  []
  [moles_H2O]
    type = ParsedFunction
    symbol_names = 'kg_H2O dt'
    symbol_values = 'kg_H2O_produced_this_timestep dt'
    expression = 'kg_H2O * 1000 / 18.0152 / dt'
  []
[]
[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 0
    bulk_modulus = 2E9
    viscosity = 1E-3
    density0 = 1000
    cv = 4000.0
    cp = 4000.0
  []
[]
[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = temperature
  mass_fraction_vars = 'f0 f1 f2'
  save_component_rate_in = 'rate_Na rate_Cl rate_SiO2 rate_H2O' # change in kg at every node / dt
  fp = the_simple_fluid
  temperature_unit = Celsius
[]
[AuxVariables]
  [injection_temperature]
    initial_condition = 200
  []
  [injection_rate_massfrac_Na]
    initial_condition = 0.002285946
  []
  [injection_rate_massfrac_Cl]
    initial_condition = 0.0035252
  []
  [injection_rate_massfrac_SiO2]
    initial_condition = 1.3741E-05
  []
  [injection_rate_massfrac_H2O]
    initial_condition = 0.994175112
  []
  [rate_H2O]
  []
  [rate_Na]
  []
  [rate_Cl]
  []
  [rate_SiO2]
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosityConst # this simulation has no porosity changes from dissolution
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0   0 1E-12 0   0 0 1E-12'
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0 0 0  0 0 0  0 0 0'
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
[]
[Preconditioning]
  active = typically_efficient
  [typically_efficient]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_hypre_type'
    petsc_options_value = ' hypre    boomeramg'
  []
  [strong]
    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      ilu           NONZERO                   2'
  []
  [probably_too_strong]
    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 = 7.76E6 # 90 days
  dt = 1E5
[]
[Outputs]
  exodus = true
[]
[MultiApps]
  [react]
    type = TransientMultiApp
    input_files = aquifer_geochemistry.i
    clone_master_mesh = true
    execute_on = 'timestep_end'
  []
[]
[Transfers]
  [changes_due_to_flow]
    type = MultiAppCopyTransfer
    source_variable = 'rate_H2O rate_Na rate_Cl rate_SiO2 temperature'
    variable = 'pf_rate_H2O pf_rate_Na pf_rate_Cl pf_rate_SiO2 temperature'
    to_multi_app = react
  []
  [massfrac_from_geochem]
    type = MultiAppCopyTransfer
    source_variable = 'massfrac_Na massfrac_Cl massfrac_SiO2'
    variable = 'f0 f1 f2'
    from_multi_app = react
  []
[]
(modules/porous_flow/test/tests/dirackernels/bh_except05.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    mass_fraction_component = 0
    point_file = bh02.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/combined/examples/geochem-porous_flow/geotes_weber_tensleep/porous_flow.i)
#########################################
#                                       #
# File written by create_input_files.py #
#                                       #
#########################################
# PorousFlow simulation of injection and production in a simplified GeoTES aquifer
# Much of this file is standard porous-flow stuff.  The unusual aspects are:
# - transfer of the rates of changes of each species (kg.s) to the aquifer_geochemistry.i simulation.  This is achieved by saving these changes from the PorousFlowMassTimeDerivative residuals
# - transfer of the temperature field to the aquifer_geochemistry.i simulation
# Interesting behaviour can be simulated by this file without its 'parent' simulation, exchanger.i.  exchanger.i provides mass-fractions injected via the injection_rate_massfrac_* variables, but since these are more-or-less constant throughout the duration of the exchanger.i simulation, the initial_conditions specified below may be used.  Similar, exchanger.i provides injection_temperature, but that is also constant.
injection_rate = -0.02 # kg/s/m, negative because injection as a source
production_rate = 0.02 # kg/s/m, this is about the maximum that can be sustained by the aquifer, with its fairly low permeability, without porepressure becoming negative
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = -75
    xmax = 75
    ymin = 0
    ymax = 40
    zmin = -25
    zmax = 25
    nx = 15
    ny = 4
    nz = 5
  []
  [aquifer]
    type = ParsedSubdomainMeshGenerator
    input = gen
    block_id = 1
    block_name = aquifer
    combinatorial_geometry = 'z >= -5 & z <= 5'
  []
  [injection_nodes]
    input = aquifer
    type = ExtraNodesetGenerator
    new_boundary = injection_nodes
    coord = '-25 0 -5; -25 0 5'
  []
  [production_nodes]
    input = injection_nodes
    type = ExtraNodesetGenerator
    new_boundary = production_nodes
    coord = '25 0 -5; 25 0 5'
  []
[]
[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 -10'
[]
[BCs]
  [injection_temperature]
    type = MatchedValueBC
    variable = temperature
    v = injection_temperature
    boundary = injection_nodes
  []
[]
[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 0
    bulk_modulus = 2E9
    viscosity = 1E-3
    density0 = 1000
    cv = 4000.0
    cp = 4000.0
  []
[]
[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = temperature
  mass_fraction_vars = 'f_H f_Cl f_SO4 f_HCO3 f_SiO2aq f_Al f_Ca f_Mg f_Fe f_K f_Na f_Sr f_F f_BOH f_Br f_Ba f_Li f_NO3 f_O2aq '
  save_component_rate_in = 'rate_H rate_Cl rate_SO4 rate_HCO3 rate_SiO2aq rate_Al rate_Ca rate_Mg rate_Fe rate_K rate_Na rate_Sr rate_F rate_BOH rate_Br rate_Ba rate_Li rate_NO3 rate_O2aq rate_H2O' # change in kg at every node / dt
  fp = the_simple_fluid
  temperature_unit = Celsius
[]
[Materials]
  [porosity_caps]
    type = PorousFlowPorosityConst # this simulation has no porosity changes from dissolution
    block = 0
    porosity = 0.01
  []
  [porosity_aquifer]
    type = PorousFlowPorosityConst # this simulation has no porosity changes from dissolution
    block = aquifer
    porosity = 0.063
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = 0
    permeability = '1E-18 0 0   0 1E-18 0   0 0 1E-18'
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1.7E-15 0 0   0 1.7E-15 0   0 0 4.1E-16'
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0 0 0  0 0 0  0 0 0'
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
[]
[Preconditioning]
  active = typically_efficient
  [typically_efficient]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_hypre_type'
    petsc_options_value = ' hypre    boomeramg'
  []
  [strong]
    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      ilu           NONZERO                   2'
  []
  [probably_too_strong]
    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 = 7.76E6 # 90 days
  [TimeStepper]
    type = FunctionDT
    function = 'min(3E4, max(1E4, 0.2 * t))'
  []
[]
[Outputs]
  exodus = true
[]
[Variables]
  [f_H]
    initial_condition = -2.952985071156e-06
  []
  [f_Cl]
    initial_condition = 0.04870664551708
  []
  [f_SO4]
    initial_condition = 0.0060359986852517
  []
  [f_HCO3]
    initial_condition = 5.0897287594019e-05
  []
  [f_SiO2aq]
    initial_condition = 3.0246609868421e-05
  []
  [f_Al]
    initial_condition = 3.268028901929e-08
  []
  [f_Ca]
    initial_condition = 0.00082159428184586
  []
  [f_Mg]
    initial_condition = 1.8546347062146e-05
  []
  [f_Fe]
    initial_condition = 4.3291908204093e-05
  []
  [f_K]
    initial_condition = 6.8434768308898e-05
  []
  [f_Na]
    initial_condition = 0.033298053919671
  []
  [f_Sr]
    initial_condition = 1.2771866652177e-05
  []
  [f_F]
    initial_condition = 5.5648860174073e-06
  []
  [f_BOH]
    initial_condition = 0.0003758574621917
  []
  [f_Br]
    initial_condition = 9.0315286107068e-05
  []
  [f_Ba]
    initial_condition = 1.5637460875161e-07
  []
  [f_Li]
    initial_condition = 8.3017067912701e-05
  []
  [f_NO3]
    initial_condition = 0.00010958455036169
  []
  [f_O2aq]
    initial_condition = -7.0806852373351e-05
  []
  [porepressure]
    initial_condition = 30E6
  []
  [temperature]
    initial_condition = 92
    scaling = 1E-6 # fluid enthalpy is roughly 1E6
  []
[]
[DiracKernels]
  [inject_H]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_H
    point_file = injection.bh
    variable = f_H
  []
  [inject_Cl]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Cl
    point_file = injection.bh
    variable = f_Cl
  []
  [inject_SO4]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_SO4
    point_file = injection.bh
    variable = f_SO4
  []
  [inject_HCO3]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_HCO3
    point_file = injection.bh
    variable = f_HCO3
  []
  [inject_SiO2aq]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_SiO2aq
    point_file = injection.bh
    variable = f_SiO2aq
  []
  [inject_Al]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Al
    point_file = injection.bh
    variable = f_Al
  []
  [inject_Ca]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Ca
    point_file = injection.bh
    variable = f_Ca
  []
  [inject_Mg]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Mg
    point_file = injection.bh
    variable = f_Mg
  []
  [inject_Fe]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Fe
    point_file = injection.bh
    variable = f_Fe
  []
  [inject_K]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_K
    point_file = injection.bh
    variable = f_K
  []
  [inject_Na]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Na
    point_file = injection.bh
    variable = f_Na
  []
  [inject_Sr]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Sr
    point_file = injection.bh
    variable = f_Sr
  []
  [inject_F]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_F
    point_file = injection.bh
    variable = f_F
  []
  [inject_BOH]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_BOH
    point_file = injection.bh
    variable = f_BOH
  []
  [inject_Br]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Br
    point_file = injection.bh
    variable = f_Br
  []
  [inject_Ba]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Ba
    point_file = injection.bh
    variable = f_Ba
  []
  [inject_Li]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_Li
    point_file = injection.bh
    variable = f_Li
  []
  [inject_NO3]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_NO3
    point_file = injection.bh
    variable = f_NO3
  []
  [inject_O2aq]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_O2aq
    point_file = injection.bh
    variable = f_O2aq
  []
  [inject_H2O]
    type = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    fluxes = ${injection_rate}
    p_or_t_vals = 0.0
    multiplying_var = injection_rate_massfrac_H2O
    point_file = injection.bh
    variable = porepressure
  []
  [produce_H]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_H
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 0
    point_file = production.bh
    variable = f_H
  []
  [produce_Cl]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Cl
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 1
    point_file = production.bh
    variable = f_Cl
  []
  [produce_SO4]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_SO4
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 2
    point_file = production.bh
    variable = f_SO4
  []
  [produce_HCO3]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_HCO3
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 3
    point_file = production.bh
    variable = f_HCO3
  []
  [produce_SiO2aq]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_SiO2aq
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 4
    point_file = production.bh
    variable = f_SiO2aq
  []
  [produce_Al]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Al
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 5
    point_file = production.bh
    variable = f_Al
  []
  [produce_Ca]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Ca
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 6
    point_file = production.bh
    variable = f_Ca
  []
  [produce_Mg]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Mg
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 7
    point_file = production.bh
    variable = f_Mg
  []
  [produce_Fe]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Fe
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 8
    point_file = production.bh
    variable = f_Fe
  []
  [produce_K]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_K
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 9
    point_file = production.bh
    variable = f_K
  []
  [produce_Na]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Na
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 10
    point_file = production.bh
    variable = f_Na
  []
  [produce_Sr]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Sr
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 11
    point_file = production.bh
    variable = f_Sr
  []
  [produce_F]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_F
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 12
    point_file = production.bh
    variable = f_F
  []
  [produce_BOH]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_BOH
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 13
    point_file = production.bh
    variable = f_BOH
  []
  [produce_Br]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Br
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 14
    point_file = production.bh
    variable = f_Br
  []
  [produce_Ba]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Ba
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 15
    point_file = production.bh
    variable = f_Ba
  []
  [produce_Li]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_Li
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 16
    point_file = production.bh
    variable = f_Li
  []
  [produce_NO3]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_NO3
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 17
    point_file = production.bh
    variable = f_NO3
  []
  [produce_O2aq]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_O2aq
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 18
    point_file = production.bh
    variable = f_O2aq
  []
  [produce_H2O]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_mass_H2O
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    mass_fraction_component = 19
    point_file = production.bh
    variable = porepressure
  []
  [produce_heat]
    type = PorousFlowPolyLineSink
    SumQuantityUO = produced_heat
    fluxes = ${production_rate}
    p_or_t_vals = 0.0
    use_enthalpy = true
    point_file = production.bh
    variable = temperature
  []
[]
[UserObjects]
  [injected_mass]
    type = PorousFlowSumQuantity
  []
  [produced_mass_H]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Cl]
    type = PorousFlowSumQuantity
  []
  [produced_mass_SO4]
    type = PorousFlowSumQuantity
  []
  [produced_mass_HCO3]
    type = PorousFlowSumQuantity
  []
  [produced_mass_SiO2aq]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Al]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Ca]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Mg]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Fe]
    type = PorousFlowSumQuantity
  []
  [produced_mass_K]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Na]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Sr]
    type = PorousFlowSumQuantity
  []
  [produced_mass_F]
    type = PorousFlowSumQuantity
  []
  [produced_mass_BOH]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Br]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Ba]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Li]
    type = PorousFlowSumQuantity
  []
  [produced_mass_NO3]
    type = PorousFlowSumQuantity
  []
  [produced_mass_O2aq]
    type = PorousFlowSumQuantity
  []
  [produced_mass_H2O]
    type = PorousFlowSumQuantity
  []
  [produced_heat]
    type = PorousFlowSumQuantity
  []
[]
[Postprocessors]
  [dt]
    type = TimestepSize
    execute_on = TIMESTEP_BEGIN
  []
  [tot_kg_injected_this_timestep]
    type = PorousFlowPlotQuantity
    uo = injected_mass
  []
  [kg_H_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_H
  []
  [kg_Cl_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Cl
  []
  [kg_SO4_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_SO4
  []
  [kg_HCO3_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_HCO3
  []
  [kg_SiO2aq_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_SiO2aq
  []
  [kg_Al_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Al
  []
  [kg_Ca_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Ca
  []
  [kg_Mg_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Mg
  []
  [kg_Fe_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Fe
  []
  [kg_K_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_K
  []
  [kg_Na_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Na
  []
  [kg_Sr_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Sr
  []
  [kg_F_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_F
  []
  [kg_BOH_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_BOH
  []
  [kg_Br_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Br
  []
  [kg_Ba_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Ba
  []
  [kg_Li_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Li
  []
  [kg_NO3_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_NO3
  []
  [kg_O2aq_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_O2aq
  []
  [kg_H2O_produced_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_mass_H2O
  []
  [mole_rate_H_produced]
    type = FunctionValuePostprocessor
    function = moles_H
    indirect_dependencies = 'kg_H_produced_this_timestep dt'
  []
  [mole_rate_Cl_produced]
    type = FunctionValuePostprocessor
    function = moles_Cl
    indirect_dependencies = 'kg_Cl_produced_this_timestep dt'
  []
  [mole_rate_SO4_produced]
    type = FunctionValuePostprocessor
    function = moles_SO4
    indirect_dependencies = 'kg_SO4_produced_this_timestep dt'
  []
  [mole_rate_HCO3_produced]
    type = FunctionValuePostprocessor
    function = moles_HCO3
    indirect_dependencies = 'kg_HCO3_produced_this_timestep dt'
  []
  [mole_rate_SiO2aq_produced]
    type = FunctionValuePostprocessor
    function = moles_SiO2aq
    indirect_dependencies = 'kg_SiO2aq_produced_this_timestep dt'
  []
  [mole_rate_Al_produced]
    type = FunctionValuePostprocessor
    function = moles_Al
    indirect_dependencies = 'kg_Al_produced_this_timestep dt'
  []
  [mole_rate_Ca_produced]
    type = FunctionValuePostprocessor
    function = moles_Ca
    indirect_dependencies = 'kg_Ca_produced_this_timestep dt'
  []
  [mole_rate_Mg_produced]
    type = FunctionValuePostprocessor
    function = moles_Mg
    indirect_dependencies = 'kg_Mg_produced_this_timestep dt'
  []
  [mole_rate_Fe_produced]
    type = FunctionValuePostprocessor
    function = moles_Fe
    indirect_dependencies = 'kg_Fe_produced_this_timestep dt'
  []
  [mole_rate_K_produced]
    type = FunctionValuePostprocessor
    function = moles_K
    indirect_dependencies = 'kg_K_produced_this_timestep dt'
  []
  [mole_rate_Na_produced]
    type = FunctionValuePostprocessor
    function = moles_Na
    indirect_dependencies = 'kg_Na_produced_this_timestep dt'
  []
  [mole_rate_Sr_produced]
    type = FunctionValuePostprocessor
    function = moles_Sr
    indirect_dependencies = 'kg_Sr_produced_this_timestep dt'
  []
  [mole_rate_F_produced]
    type = FunctionValuePostprocessor
    function = moles_F
    indirect_dependencies = 'kg_F_produced_this_timestep dt'
  []
  [mole_rate_BOH_produced]
    type = FunctionValuePostprocessor
    function = moles_BOH
    indirect_dependencies = 'kg_BOH_produced_this_timestep dt'
  []
  [mole_rate_Br_produced]
    type = FunctionValuePostprocessor
    function = moles_Br
    indirect_dependencies = 'kg_Br_produced_this_timestep dt'
  []
  [mole_rate_Ba_produced]
    type = FunctionValuePostprocessor
    function = moles_Ba
    indirect_dependencies = 'kg_Ba_produced_this_timestep dt'
  []
  [mole_rate_Li_produced]
    type = FunctionValuePostprocessor
    function = moles_Li
    indirect_dependencies = 'kg_Li_produced_this_timestep dt'
  []
  [mole_rate_NO3_produced]
    type = FunctionValuePostprocessor
    function = moles_NO3
    indirect_dependencies = 'kg_NO3_produced_this_timestep dt'
  []
  [mole_rate_O2aq_produced]
    type = FunctionValuePostprocessor
    function = moles_O2aq
    indirect_dependencies = 'kg_O2aq_produced_this_timestep dt'
  []
  [mole_rate_H2O_produced]
    type = FunctionValuePostprocessor
    function = moles_H2O
    indirect_dependencies = 'kg_H2O_produced_this_timestep dt'
  []
  [heat_joules_extracted_this_timestep]
    type = PorousFlowPlotQuantity
    uo = produced_heat
  []
  [production_temperature]
    type = AverageNodalVariableValue
    boundary = production_nodes
    variable = temperature
  []
[]
[Functions]
  [moles_H]
    type = ParsedFunction
    symbol_names = 'kg_H dt'
    symbol_values = 'kg_H_produced_this_timestep dt'
    expression = 'kg_H * 1000 / 1.0079 / dt'
  []
  [moles_Cl]
    type = ParsedFunction
    symbol_names = 'kg_Cl dt'
    symbol_values = 'kg_Cl_produced_this_timestep dt'
    expression = 'kg_Cl * 1000 / 35.453 / dt'
  []
  [moles_SO4]
    type = ParsedFunction
    symbol_names = 'kg_SO4 dt'
    symbol_values = 'kg_SO4_produced_this_timestep dt'
    expression = 'kg_SO4 * 1000 / 96.0576 / dt'
  []
  [moles_HCO3]
    type = ParsedFunction
    symbol_names = 'kg_HCO3 dt'
    symbol_values = 'kg_HCO3_produced_this_timestep dt'
    expression = 'kg_HCO3 * 1000 / 61.0171 / dt'
  []
  [moles_SiO2aq]
    type = ParsedFunction
    symbol_names = 'kg_SiO2aq dt'
    symbol_values = 'kg_SiO2aq_produced_this_timestep dt'
    expression = 'kg_SiO2aq * 1000 / 60.0843 / dt'
  []
  [moles_Al]
    type = ParsedFunction
    symbol_names = 'kg_Al dt'
    symbol_values = 'kg_Al_produced_this_timestep dt'
    expression = 'kg_Al * 1000 / 26.9815 / dt'
  []
  [moles_Ca]
    type = ParsedFunction
    symbol_names = 'kg_Ca dt'
    symbol_values = 'kg_Ca_produced_this_timestep dt'
    expression = 'kg_Ca * 1000 / 40.08 / dt'
  []
  [moles_Mg]
    type = ParsedFunction
    symbol_names = 'kg_Mg dt'
    symbol_values = 'kg_Mg_produced_this_timestep dt'
    expression = 'kg_Mg * 1000 / 24.305 / dt'
  []
  [moles_Fe]
    type = ParsedFunction
    symbol_names = 'kg_Fe dt'
    symbol_values = 'kg_Fe_produced_this_timestep dt'
    expression = 'kg_Fe * 1000 / 55.847 / dt'
  []
  [moles_K]
    type = ParsedFunction
    symbol_names = 'kg_K dt'
    symbol_values = 'kg_K_produced_this_timestep dt'
    expression = 'kg_K * 1000 / 39.0983 / dt'
  []
  [moles_Na]
    type = ParsedFunction
    symbol_names = 'kg_Na dt'
    symbol_values = 'kg_Na_produced_this_timestep dt'
    expression = 'kg_Na * 1000 / 22.9898 / dt'
  []
  [moles_Sr]
    type = ParsedFunction
    symbol_names = 'kg_Sr dt'
    symbol_values = 'kg_Sr_produced_this_timestep dt'
    expression = 'kg_Sr * 1000 / 87.62 / dt'
  []
  [moles_F]
    type = ParsedFunction
    symbol_names = 'kg_F dt'
    symbol_values = 'kg_F_produced_this_timestep dt'
    expression = 'kg_F * 1000 / 18.9984 / dt'
  []
  [moles_BOH]
    type = ParsedFunction
    symbol_names = 'kg_BOH dt'
    symbol_values = 'kg_BOH_produced_this_timestep dt'
    expression = 'kg_BOH * 1000 / 61.8329 / dt'
  []
  [moles_Br]
    type = ParsedFunction
    symbol_names = 'kg_Br dt'
    symbol_values = 'kg_Br_produced_this_timestep dt'
    expression = 'kg_Br * 1000 / 79.904 / dt'
  []
  [moles_Ba]
    type = ParsedFunction
    symbol_names = 'kg_Ba dt'
    symbol_values = 'kg_Ba_produced_this_timestep dt'
    expression = 'kg_Ba * 1000 / 137.33 / dt'
  []
  [moles_Li]
    type = ParsedFunction
    symbol_names = 'kg_Li dt'
    symbol_values = 'kg_Li_produced_this_timestep dt'
    expression = 'kg_Li * 1000 / 6.941 / dt'
  []
  [moles_NO3]
    type = ParsedFunction
    symbol_names = 'kg_NO3 dt'
    symbol_values = 'kg_NO3_produced_this_timestep dt'
    expression = 'kg_NO3 * 1000 / 62.0049 / dt'
  []
  [moles_O2aq]
    type = ParsedFunction
    symbol_names = 'kg_O2aq dt'
    symbol_values = 'kg_O2aq_produced_this_timestep dt'
    expression = 'kg_O2aq * 1000 / 31.9988 / dt'
  []
  [moles_H2O]
    type = ParsedFunction
    symbol_names = 'kg_H2O dt'
    symbol_values = 'kg_H2O_produced_this_timestep dt'
    expression = 'kg_H2O * 1000 / 18.01801802 / dt'
  []
[]
[AuxVariables]
  [injection_temperature]
    initial_condition = 92
  []
  [injection_rate_massfrac_H]
    initial_condition = -2.952985071156e-06
  []
  [injection_rate_massfrac_Cl]
    initial_condition = 0.04870664551708
  []
  [injection_rate_massfrac_SO4]
    initial_condition = 0.0060359986852517
  []
  [injection_rate_massfrac_HCO3]
    initial_condition = 5.0897287594019e-05
  []
  [injection_rate_massfrac_SiO2aq]
    initial_condition = 3.0246609868421e-05
  []
  [injection_rate_massfrac_Al]
    initial_condition = 3.268028901929e-08
  []
  [injection_rate_massfrac_Ca]
    initial_condition = 0.00082159428184586
  []
  [injection_rate_massfrac_Mg]
    initial_condition = 1.8546347062146e-05
  []
  [injection_rate_massfrac_Fe]
    initial_condition = 4.3291908204093e-05
  []
  [injection_rate_massfrac_K]
    initial_condition = 6.8434768308898e-05
  []
  [injection_rate_massfrac_Na]
    initial_condition = 0.033298053919671
  []
  [injection_rate_massfrac_Sr]
    initial_condition = 1.2771866652177e-05
  []
  [injection_rate_massfrac_F]
    initial_condition = 5.5648860174073e-06
  []
  [injection_rate_massfrac_BOH]
    initial_condition = 0.0003758574621917
  []
  [injection_rate_massfrac_Br]
    initial_condition = 9.0315286107068e-05
  []
  [injection_rate_massfrac_Ba]
    initial_condition = 1.5637460875161e-07
  []
  [injection_rate_massfrac_Li]
    initial_condition = 8.3017067912701e-05
  []
  [injection_rate_massfrac_NO3]
    initial_condition = 0.00010958455036169
  []
  [injection_rate_massfrac_O2aq]
    initial_condition = -7.0806852373351e-05
  []
  [injection_rate_massfrac_H2O]
    initial_condition = 0.91032275033842
  []
  [rate_H]
  []
  [rate_Cl]
  []
  [rate_SO4]
  []
  [rate_HCO3]
  []
  [rate_SiO2aq]
  []
  [rate_Al]
  []
  [rate_Ca]
  []
  [rate_Mg]
  []
  [rate_Fe]
  []
  [rate_K]
  []
  [rate_Na]
  []
  [rate_Sr]
  []
  [rate_F]
  []
  [rate_BOH]
  []
  [rate_Br]
  []
  [rate_Ba]
  []
  [rate_Li]
  []
  [rate_NO3]
  []
  [rate_O2aq]
  []
  [rate_H2O]
  []
[]
[MultiApps]
  [react]
    type = TransientMultiApp
    input_files = aquifer_geochemistry.i
    clone_master_mesh = true
    execute_on = 'timestep_end'
  []
[]
[Transfers]
  [changes_due_to_flow]
    type = MultiAppCopyTransfer
    source_variable = 'rate_H rate_Cl rate_SO4 rate_HCO3 rate_SiO2aq rate_Al rate_Ca rate_Mg rate_Fe rate_K rate_Na rate_Sr rate_F rate_BOH rate_Br rate_Ba rate_Li rate_NO3 rate_O2aq rate_H2O temperature'
    variable = 'pf_rate_H pf_rate_Cl pf_rate_SO4 pf_rate_HCO3 pf_rate_SiO2aq pf_rate_Al pf_rate_Ca pf_rate_Mg pf_rate_Fe pf_rate_K pf_rate_Na pf_rate_Sr pf_rate_F pf_rate_BOH pf_rate_Br pf_rate_Ba pf_rate_Li pf_rate_NO3 pf_rate_O2aq pf_rate_H2O temperature'
    to_multi_app = react
  []
  [massfrac_from_geochem]
    type = MultiAppCopyTransfer
    source_variable = 'massfrac_H massfrac_Cl massfrac_SO4 massfrac_HCO3 massfrac_SiO2aq massfrac_Al massfrac_Ca massfrac_Mg massfrac_Fe massfrac_K massfrac_Na massfrac_Sr massfrac_F massfrac_BOH massfrac_Br massfrac_Ba massfrac_Li massfrac_NO3 massfrac_O2aq '
    variable = 'f_H f_Cl f_SO4 f_HCO3 f_SiO2aq f_Al f_Ca f_Mg f_Fe f_K f_Na f_Sr f_F f_BOH f_Br f_Ba f_Li f_NO3 f_O2aq '
    from_multi_app = react
  []
[]
(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
  []
[]
[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/dirackernels/bh_except06.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    use_mobility = true
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/pls01.i)
# fully-saturated situation with a poly-line sink at one
# of the nodes.  Because there is no fluid flow, the
# other nodes should not experience any change in
# porepressure.
# The poly-line sink has length=2 and weight=0.1, and
# extracts fluid at a constant rate of 1 kg.m^-1.s^-1.
# Therefore, in 1 second it will have extracted a total
# of 0.2 kg.
# The porosity is 0.1, and the elemental volume is 2,
# so the fluid mass at the node in question = 0.2 * density / 4,
# where the 4 is the number of nodes in the element.
# In this simulation density = dens0 * exp(P / bulk), with
# dens0 = 100, and bulk = 20 MPa.
# The initial porepressure P0 = 10 MPa, so the final (after
# 1 second of simulation) is
# P(t=1) = 0.950879 MPa
#
[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 0
  xmax = 2
  ymin = 0
  ymax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    density0 = 100
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]
[DiracKernels]
  [pls]
    type = PorousFlowPolyLineSink
    fluid_phase = 0
    point_file = pls01_21.bh
    line_length = 2
    SumQuantityUO = pls_total_outflow_mass
    variable = pp
    p_or_t_vals = '0 1E7'
    fluxes = '1 1'
  []
[]
[Postprocessors]
  [pls_report]
    type = PorousFlowPlotQuantity
    uo = pls_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 pls_report'
  []
  [p00]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
  [p01]
    type = PointValue
    variable = pp
    point = '0 1 0'
    execute_on = timestep_end
  []
  [p20]
    type = PointValue
    variable = pp
    point = '2 0 0'
    execute_on = timestep_end
  []
  [p21]
    type = PointValue
    variable = pp
    point = '2 1 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 pls_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]
[Outputs]
  file_base = pls01
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except01.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 1
    point_file = bh02.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh04.i)
# fully-saturated
# production
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Functions]
  [dts]
    type = PiecewiseLinear
    y = '1E-2 1E-1 1 1E1 1E2 1E3'
    x = '0 1E-1 1 1E1 1E2 1E3'
  []
[]
[Variables]
  [pp]
    initial_condition = 0
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.8
    alpha = 1e-5
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC
    m = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh02.bh
    fluid_phase = 0
    bottom_p_or_t = -1E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 bh_report'
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 1E3
  solve_type = NEWTON
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]
[Outputs]
  file_base = bh04
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except14.i)
# fully-saturated
# production
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02_huge.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/actions/basicthm_borehole.i)
# PorousFlowBasicTHM action with coupling_type = Hydro (no thermal or
# mechanical effects), plus a Peaceman borehole with use_mobility = true
# to test that nodal relative permeability is added by this action.
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [porepressure]
    initial_condition = 1e7
  []
[]
[AuxVariables]
  [temperature]
    initial_condition = 293
  []
[]
[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = simple_fluid
  multiply_by_density = true
[]
[UserObjects]
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 1
    fluid_bulk_modulus = 2e9
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0   0 1e-13 0   0 0 1e-13'
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = borehole_total_outflow_mass
    point_file = borehole.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 0
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1e-10 1e-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 0.1
  solve_type = NEWTON
[]
[Outputs]
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/actions/fullsat_borehole.i)
# PorousFlowFullySaturated action with coupling_type = ThermoHydro (no
# mechanical effects), plus a Peaceman borehole with use_mobility = true
# to test that nodal relative permeability is added by this action.
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [porepressure]
    initial_condition = 1E7
  []
  [temperature]
    initial_condition = 323.15
  []
[]
[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = temperature
  dictator_name = dictator
  stabilization = none
  fp = simple_fluid
  gravity = '0 0 0'
[]
[BCs]
  [temperature]
    type = DirichletBC
    variable = temperature
    boundary = 'left right'
    value = 323.15
  []
[]
[UserObjects]
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
  []
[]
[Materials]
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.25
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0   0 1e-13 0   0 0 1e-13'
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '3 0 0  0 3 0  0 0 3'
    wet_thermal_conductivity = '3 0 0  0 3 0  0 0 3'
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 850
    density = 2700
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = borehole_total_outflow_mass
    point_file = borehole.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 0
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 0.1
  solve_type = NEWTON
[]
[Outputs]
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except12.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = does_not_exist
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh05.i)
# unsaturated
# injection
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Functions]
  [dts]
    type = PiecewiseLinear
    y = '500 500 1E1'
    x = '4000 5000 6500'
  []
[]
[Variables]
  [pp]
    initial_condition = -2E5
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.8
    alpha = 1e-5
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC
    m = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh03.bh
    fluid_phase = 0
    bottom_p_or_t = 0
    unit_weight = '0 0 0'
    use_mobility = true
    character = -1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 bh_report'
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 6500
  solve_type = NEWTON
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
[]
[Outputs]
  file_base = bh05
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh03.i)
# fully-saturated
# injection
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 0
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    variable = pp
    SumQuantityUO = borehole_total_outflow_mass
    point_file = bh03.bh
    function_of = pressure
    fluid_phase = 0
    bottom_p_or_t = 'insitu_pp'
    unit_weight = '0 0 0'
    use_mobility = true
    character = -1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 bh_report'
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
  [insitu_pp]
    type = ParsedFunction
    expression = '-2e7*z' #bh_bottom is located at z=-0.5
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
[Outputs]
  file_base = bh03
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh02reporter.i)
# fully-saturated
# production
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    # Because the Variable for this Sink is pp, and pp is associated
    # with the fluid-mass conservation equation, this sink is extracting
    # fluid mass (and not heat energy or something else)
    variable = pp
    # The following specfies that the total fluid mass coming out of
    # the porespace via this sink in this timestep should be recorded
    # in the pls_total_outflow_mass UserObject
    SumQuantityUO = borehole_total_outflow_mass
    # The following file defines the polyline geometry
    # which is just two points in this particular example
    weight_reporter='bh02file/column_0'
    x_coord_reporter='bh02file/column_1'
    y_coord_reporter='bh02file/column_2'
    z_coord_reporter='bh02file/column_3'
    # First, we want Peacemans f to be a function of porepressure (and not
    # temperature or something else).  So bottom_p_or_t is actually porepressure
    function_of = pressure
    fluid_phase = 0
    # The bottomhole pressure
    bottom_p_or_t = 0
    # In this example there is no increase of the wellbore pressure
    # due to gravity:
    unit_weight = '0 0 0'
    # PeacemanBoreholes should almost always have use_mobility = true
    use_mobility = true
    # This is a production wellbore (a sink of fluid that removes fluid from porespace)
    character = 1
  []
[]
[VectorPostprocessors]
  [bh02file]
    type = CSVReader
    csv_file = bh02.bh
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
    indirect_dependencies = 'fluid_mass1 fluid_mass0 bh_report'
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
[Outputs]
  exodus = false
  csv = true
  execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except15.i)
# fully-saturated
# production
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/injection_with_plasticity.i)
# Example: Injection into a uniform aquifer 10 x 10 x 5 km
# Drucker-Prager deformation
# Darcy flow
gravity = -9.81
solid_density = 2350
fluid_density = 1000
porosity0 = 0.1
[Mesh]
  type = GeneratedMesh
  dim = 3
  xmin = 0
  xmax = 1e4
  ymin = 0
  ymax = 1e4
  zmax = 0
  zmin = -5e3
  nx = 2
  ny = 2
  nz = 2
[]
[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 ${gravity}'
  displacements = 'disp_x disp_y disp_z'
  strain_at_nearest_qp = true
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 0 # Not doing a thermal simulation
    bulk_modulus = 2E9
    density0 = ${fluid_density}
    viscosity = 5E-4
  []
[]
[PorousFlowFullySaturated]
  coupling_type = HydroMechanical
  porepressure = pp
  dictator_name = dictator
  fp = simple_fluid
  add_darcy_aux = false
  add_stress_aux = false
  stabilization = none
[]
[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [pp]
    scaling = 1E6
    [InitialCondition]
      type = FunctionIC
      function = ini_pp
    []
  []
[]
[Functions]
  [ini_stress]
    type = ParsedFunction
    expression = '-${gravity} * z * (${solid_density} - ${fluid_density}) * (1.0 - ${porosity0})'  # initial effective stress that should result from weight force
  []
  [ini_pp]
    type = ParsedFunction
    expression = '${gravity} * z * ${fluid_density} + 1E5'
  []
[]
[BCs]
  [p_top]
    type = FunctionDirichletBC
    variable = pp
    boundary = front
    function = ini_pp
  []
  [x_roller]
    type = DirichletBC
    variable = disp_x
    boundary = 'left right'
    value = 0
  []
  [y_roller]
    type = DirichletBC
    variable = disp_y
    boundary = 'top bottom'
    value = 0
  []
  [z_confined]
    type = DirichletBC
    variable = disp_z
    boundary = 'back front'
    value = 0
  []
[]
[UserObjects]
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  # Cohesion
  [mc_coh]
    type = TensorMechanicsHardeningConstant
    value = 6.0E6
  []
  # Friction angle
  [mc_phi]
    type = TensorMechanicsHardeningConstant
    value = 35.0
    convert_to_radians = true
  []
  # Dilation angle
  [mc_psi]
    type = TensorMechanicsHardeningConstant
    value = 2
    convert_to_radians = true
  []
  # Drucker-Prager objects
  [dp]
    type = TensorMechanicsPlasticDruckerPragerHyperbolic
    mc_cohesion = mc_coh
    mc_friction_angle = mc_phi
    mc_dilation_angle = mc_psi
    yield_function_tolerance = 1E-3
    internal_constraint_tolerance = 1E-6
  []
  # Tensile strength
  [tens]
    type = TensorMechanicsHardeningConstant
    value = 3.0E6
  []
  # Compressive strength (cap on yield envelope)
  [compr_all]
    type = TensorMechanicsHardeningConstant
    value = 1E10
  []
[]
[Materials]
  [strain]
    type = ComputeIncrementalStrain
    eigenstrain_names = eigenstrain_all
  []
  [eigenstrain_all]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = 'ini_stress 0 0  0 ini_stress 0  0 0 ini_stress'
    eigenstrain_name = eigenstrain_all
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    bulk_modulus = 3.3333E9
    shear_modulus = 2.5E9
  []
  [dp_mat]
    type = CappedDruckerPragerStressUpdate
    DP_model = dp
    tensile_strength = tens
    compressive_strength = compr_all
    smoothing_tol = 1E5
    yield_function_tol = 1E-3
    tip_smoother = 0
  []
  [stress]
    type = ComputeMultipleInelasticStress
    inelastic_models = dp_mat
  []
  # Permeability
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-13 0 0  0 1E-13 0  0 0 1E-13'
  []
  # Porosity
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = ${porosity0}
    biot_coefficient = 1.0
    solid_bulk = 1.0 # Required but irrelevant when biot_coefficient is unity
    mechanical = true
    fluid = true
  []
  # Density of saturated rock
  [density]
    type = PorousFlowTotalGravitationalDensityFullySaturatedFromPorosity
    rho_s = ${solid_density}
  []
[]
[DiracKernels]
  [pls]
    type = PorousFlowPolyLineSink
    variable = pp
    SumQuantityUO = pls_total_outflow_mass
    point_file = two_nodes.bh
    function_of = pressure
    fluid_phase = 0
    p_or_t_vals = '0 1E7'
    fluxes = '-1.59 -1.59'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
  []
[]
[Executioner]
  solve_type = Newton
  type = Transient
  dt = 1E6
  end_time = 1E6
  nl_rel_tol = 1E-7
[]
[Outputs]
  exodus = true
[]
(modules/porous_flow/test/tests/dirackernels/bh_except07.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]
[GlobalParams]
  PorousFlowDictator = dictator
[]
[Variables]
  [pp]
    initial_condition = 1E7
  []
[]
[Kernels]
  [mass0]
    type = TimeDerivative
    variable = pp
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]
[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]
[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    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 = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
[]
[DiracKernels]
  [bh]
    type = PorousFlowPeacemanBorehole
    bottom_p_or_t = 0
    fluid_phase = 0
    point_file = bh02.bh
    use_mobility = true
    SumQuantityUO = borehole_total_outflow_mass
    variable = pp
    unit_weight = '0 0 0'
    character = 1
  []
[]
[Postprocessors]
  [bh_report]
    type = PorousFlowPlotQuantity
    uo = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass
    execute_on = timestep_begin
  []
  [fluid_mass1]
    type = PorousFlowFluidMass
    execute_on = timestep_end
  []
  [zmass_error]
    type = FunctionValuePostprocessor
    function = mass_bal_fcn
    execute_on = timestep_end
  []
  [p0]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = timestep_end
  []
[]
[Functions]
  [mass_bal_fcn]
    type = ParsedFunction
    expression = abs((a-c+d)/2/(a+c))
    symbol_names = 'a c d'
    symbol_values = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]
[Preconditioning]
  [usual]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]
[Executioner]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]