Point and line sources and sinks

This page may be read in conjunction with the theory behind Dirac Kernels in PorousFlow.

Geometric tests

The test suite contains many core tests that demonstrate:

when a point sink is placed at a node, it withdraws fluid (or heat) from only that node, for instance

# 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<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 0
  xmax = 2
  ymin = 0
  ymax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1E7
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../../../source/userobjects/PorousFlowSumQuantity.html"}>>>
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.5
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-7
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e7
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 100
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.1
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [pls]
    type = PorousFlowPolyLineSink<<<{"description": "Approximates a polyline sink by using a number of point sinks with given weighting whose positions are read from a file.  NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../../../source/dirackernels/PorousFlowPolyLineSink.html"}>>>
    fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink.  For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0
    point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink.  The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii.  Each line in the file must contain a space-separated weight and coordinate, viz r x y z.  For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure.  If your file contains just one point, you must also specify the line_length and line_direction parameters.  Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = pls01_21.bh
    line_length<<<{"description": "Line length.  Note this is only used if there is only one point in the point_file."}>>> = 2
    SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = pls_total_outflow_mass
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    p_or_t_vals<<<{"description": "Tuple of pressure (or temperature) values.  Must be monotonically increasing."}>>> = '0 1E7'
    fluxes<<<{"description": "Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed).  These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1.  A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature).  If a quad-point pressure is less than the first pressure value, the first flux value is used.  If quad-point pressure exceeds the final pressure value, the final flux value is used.  This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE."}>>> = '1 1'
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [pls_report]
    type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
    uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = pls_total_outflow_mass
  []

  [fluid_mass0]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_begin
  []

  [fluid_mass1]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_end
  []

  [zmass_error]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = mass_bal_fcn
    execute_on<<<{"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."}>>> = timestep_end
    indirect_dependencies<<<{"description": "If the evaluated function depends on other postprocessors they must be listed here to ensure proper dependency resolution"}>>> = 'fluid_mass1 fluid_mass0 pls_report'
  []

  [p00]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
  [p01]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 1 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
  [p20]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '2 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
  [p21]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '2 1 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [mass_bal_fcn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = abs((a-c+d)/2/(a+c))
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'a c d'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'fluid_mass1 fluid_mass0 pls_report'
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [usual]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = pls01
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"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."}>>> = timestep_end
[]

when a point sink that is proportional to mobility (or relative permeability, etc) is placed in an element where some nodes have zero mobility (or relative permeability, etc), then fluid (or heat) is not extracted from those nodes, for instance

# 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<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 1
  ny = 1
  xmin = 0
  xmax = 2
  ymin = 0
  ymax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
  []
[]

[ICs<<<{"href": "../../../../syntax/ICs/index.html"}>>>]
  [pp]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pp
    #function = if((x<1)&(y<0.5),1E7,-1E7)
    function<<<{"description": "The initial condition 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<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pls_total_outflow_mass]
    type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../../../source/userobjects/PorousFlowSumQuantity.html"}>>>
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.5
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-7
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e7
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 100
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.1
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC<<<{"description": "This Material calculates relative permeability of a phase using a model inspired by FLAC", "href": "../../../../source/materials/PorousFlowRelativePermeabilityFLAC.html"}>>>
    phase<<<{"description": "The phase number"}>>> = 0
    m<<<{"description": "relperm = (1 + m)seff^m - m seff^(m+1)"}>>> = 2
    s_res<<<{"description": "The residual saturation of the phase j. Must be between 0 and 1"}>>> = 0.99
    sum_s_res<<<{"description": "Sum of residual saturations over all phases.  Must be between 0 and 1"}>>> = 0.99
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [pls]
    type = PorousFlowPolyLineSink<<<{"description": "Approximates a polyline sink by using a number of point sinks with given weighting whose positions are read from a file.  NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../../../source/dirackernels/PorousFlowPolyLineSink.html"}>>>
    fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink.  For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0
    point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink.  The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii.  Each line in the file must contain a space-separated weight and coordinate, viz r x y z.  For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure.  If your file contains just one point, you must also specify the line_length and line_direction parameters.  Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = pls03.bh
    use_relative_permeability<<<{"description": "Multiply the flux by the fluid relative permeability"}>>> = true
    line_length<<<{"description": "Line length.  Note this is only used if there is only one point in the point_file."}>>> = 4
    SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = pls_total_outflow_mass
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    p_or_t_vals<<<{"description": "Tuple of pressure (or temperature) values.  Must be monotonically increasing."}>>> = '0 1E7'
    fluxes<<<{"description": "Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed).  These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1.  A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature).  If a quad-point pressure is less than the first pressure value, the first flux value is used.  If quad-point pressure exceeds the final pressure value, the final flux value is used.  This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE."}>>> = '1 1'
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [pls_report]
    type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
    uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = pls_total_outflow_mass
  []

  [fluid_mass0]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_begin
  []

  [fluid_mass1]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_end
  []

  [zmass_error]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = mass_bal_fcn
    execute_on<<<{"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."}>>> = timestep_end
    indirect_dependencies<<<{"description": "If the evaluated function depends on other postprocessors they must be listed here to ensure proper dependency resolution"}>>> = 'fluid_mass1 fluid_mass0 pls_report'
  []

  [p00]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
  [p01]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 1 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
  [p20]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '2 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
  [p21]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '2 1 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [mass_bal_fcn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = abs((a-c+d)/2/(a+c))
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'a c d'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'fluid_mass1 fluid_mass0 pls_report'
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [usual]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  end_time = 1
  dt = 1
  solve_type = NEWTON
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = pls03
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"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."}>>> = timestep_end
[]

Basic point sources/sinks

The following input file extracts fluid at the rate (kg.s $var element = document.getElementById("moose-equation-2a7dd690-1b22-4799-9dd7-d904e4a0b728");katex.render("^{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ ) $var element = document.getElementById("moose-equation-bd1259c3-9bfc-4881-add3-bb4807abdfdf");katex.render("\\mathrm{source} = \\left\\{ \\begin{array}{ll} -0.1 & \\mathrm{\\ for\\ } 100 \\leq t \\leq 300 \\\\ -0.1 & \\mathrm{\\ for\\ } 600 \\leq t \\leq 1400 \\\\ 0.2 & \\mathrm{\\ for\\ } 1500 \\leq t \\leq 2000 \\\\ 0 & \\mathrm{\\ otherwise} \\end{array} \\right.", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

# Test PorousFlowSquarePulsePointSource DiracKernel

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  bias_x = 1.1
  bias_y = 1.1
  ymax = 1
  xmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = pp
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.5
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-7
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.2
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [total_mass]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = 'initial timestep_end'
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  nl_abs_tol = 1e-14
  dt = 200
  end_time = 2000
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = squarepulse1
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"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."}>>> = 'initial timestep_end'
  [con]
    output_linear<<<{"description": "Specifies whether output occurs on each PETSc linear residual evaluation"}>>> = true
    type = Console<<<{"description": "Object for screen output.", "href": "../../../../source/outputs/Console.html"}>>>
  []
[]

[ICs<<<{"href": "../../../../syntax/ICs/index.html"}>>>]
  [PressureIC]
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pp
    type = ConstantIC<<<{"description": "Sets a constant field value.", "href": "../../../../source/ics/ConstantIC.html"}>>>
    value<<<{"description": "The value to be set in IC"}>>> = 20e6
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [sink1]
    type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
    start_time<<<{"description": "The time at which the source will start (Default is 0)"}>>> = 100
    end_time<<<{"description": "The time at which the source will end (Default is 1e30)"}>>> = 300
    point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0.5 0.5 0'
    mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = -0.1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
  [sink]
    type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
    start_time<<<{"description": "The time at which the source will start (Default is 0)"}>>> = 600
    end_time<<<{"description": "The time at which the source will end (Default is 1e30)"}>>> = 1400
    point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0.5 0.5 0'
    mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = -0.1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
  [source]
    point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0.5 0.5 0'
    start_time<<<{"description": "The time at which the source will start (Default is 0)"}>>> = 1500
    mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = 0.2
    end_time<<<{"description": "The time at which the source will end (Default is 1e30)"}>>> = 2000
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
  []
[]

MOOSE produces the expected result:

Figure 1: Results of a "square-pulsed" extraction and injection of fluid from a fluid-flow simulation.

Theis tests

In the fully-saturated, isothermal case with no mechanical coupling and constant, large fluid bulk modulus, the fluid flow equations reduce to the form conventionally used in groundwater flow: $var element = document.getElementById("moose-equation-390b6bb2-c145-4f1a-b397-f36e3f538a1a");katex.render("S_{s}\\frac{\\partial H}{\\partial t} = \\nabla_{i}\\left(K_{ij}\\nabla_{j} H \\right) \\ ,", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ where

$var element = document.getElementById("moose-equation-d2e1a0c0-e5d0-4467-91b7-ec28ed0b2155");katex.render("H", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the "head", $var element = document.getElementById("moose-equation-3b99298c-4108-4667-959c-0c6b8e1b66a8");katex.render("H=(P - P_{0})/(\\rho g) - \\mathrm{depth}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . Here $var element = document.getElementById("moose-equation-701a47f0-6e54-40cc-8908-2c425193bbd1");katex.render("\\rho", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the fluid density.
$var element = document.getElementById("moose-equation-0b8a387d-21a4-4e90-820c-ca147ca2878e");katex.render("S_{s}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the "specific storage". If the rock compressibility is ignored then $var element = document.getElementById("moose-equation-c9dc6b5f-d49e-45ae-a745-c5fae9a91e4a");katex.render("S_{s} = \\rho g \\phi / K_{\\mathrm{f}}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , where $var element = document.getElementById("moose-equation-c3db7f60-c66e-4d73-95bd-e838fe5715b9");katex.render("K_{\\mathrm{f}}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the fluid bulk modulus
$var element = document.getElementById("moose-equation-33235e2e-c5f6-469c-8e81-76acc3ab1a19");katex.render("K_{ij}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the hydraulic conductivity: $var element = document.getElementById("moose-equation-9fb55143-e163-48f7-82c3-3a2848b55850");katex.render("K_{ij} = k_{ij}\\rho g / \\mu", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , where $var element = document.getElementById("moose-equation-ea08372e-b491-4522-ab9a-d6c6a9dc5931");katex.render("k_{ij}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the permeability and $var element = document.getElementById("moose-equation-d541f880-0930-4e6b-9280-8fe08dff0310");katex.render("\\mu", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the fluid viscosity

This equation is identical to the physics described by the PorousFlowFullySaturatedMassTimeDerivative and the PorousFlowFullySaturatedDarcyBase Kernels, which may be used explicitly, or through the PorousFlowBasicTHM Action. When using these, the specific storage $var element = document.getElementById("moose-equation-41783b5f-4baa-4999-87b2-b920dfb34668");katex.render("S_{s}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ in the above equation is replaced by $var element = document.getElementById("moose-equation-81afee6f-b2e6-4592-854e-49aa931b4da3");katex.render("1/M", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ (reciprocal of the Biot modulus), and $var element = document.getElementById("moose-equation-bd98bd7f-caaf-489e-832c-8f22fe14a994");katex.render("K_{ij}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ by $var element = document.getElementById("moose-equation-9d2431d8-08e9-4594-b685-291e292d7deb");katex.render("k_{ij}/\\mu", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , and $var element = document.getElementById("moose-equation-f68bbc8b-4623-4240-9225-396979a41d72");katex.render("H", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ by $var element = document.getElementById("moose-equation-7061c6aa-a5b6-4697-9634-85dfdc0dbf3b");katex.render("P", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . Finally, remember that this formulation is volume based, not mass based like the remainder of PorousFlow.

Place a constant volumetric sink of strength $var element = document.getElementById("moose-equation-d4a1ea74-9eac-49ec-bee4-2ae1d47bce5e");katex.render("q", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ (m $var element = document.getElementById("moose-equation-55d4592d-7be3-446c-ba48-cc21cd9b8c94");katex.render("^{3}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .m $var element = document.getElementById("moose-equation-ed9742a6-b13b-431e-9ac3-046ee4e31059");katex.render("^{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .s $var element = document.getElementById("moose-equation-263957cf-8b32-495e-95e8-adde62321074");katex.render("^{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ ) acting along an infinite line in an isotropic 3D medium. The situation is therefore two dimensional. Theis provided the solution for the head $(1)var element = document.getElementById("moose-equation-db85b848-d2fd-4949-ad23-5bfe69cc051f");katex.render("H = H_{0} + \\frac{q}{4\\pi K}\\mathrm{Ei}\\left( -\\frac{r^{2}S_{s}}{4Kt} \\right) \\ , ", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ which is frequently used in the groundwater literature. "Ei" is the exponential integral function, with values $var element = document.getElementById("moose-equation-577340b6-c087-44e2-b12a-c918a5fc89a3");katex.render("\\mathrm{Ei}(-x) = \\gamma + log(x) - x + O(x^{2})", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ for small $var element = document.getElementById("moose-equation-0a041258-7240-4fe3-a890-15743cd9912a");katex.render("x", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ (where $var element = document.getElementById("moose-equation-d0dd1f6e-60d1-4379-bdd3-c4b6577b85b7");katex.render("\\gamma", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the Euler number $var element = document.getElementById("moose-equation-e0a4fb74-1d9c-4bbc-82e4-5f7ccb7f5887");katex.render("0.57722\\ldots", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ ), and $var element = document.getElementById("moose-equation-8d9da955-b80e-4ebb-81f9-e0dd63d51a2c");katex.render("\\mathrm{Ei}(-x) = e^{-x}\\left(-\\frac{1}{x} + \\frac{1}{x^{2}} + O(\\frac{1}{x^{3}})\\right)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ for large $var element = document.getElementById("moose-equation-f7bd5bb1-bcd9-43a6-b826-0bc842e02dbc");katex.render("x", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .

This is checked using a number of PorousFlow tests (the extraction specified using mass rate, with volume rate, and from single and 2-phase systems). For instance:

# Theis problem: Flow to single sink using BasicTHM
# SinglePhase
# RZ mesh

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 1
  nx = 20
  xmax = 100
  bias_x = 1.05
  coord_type = RZ
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20E6
  []
[]

[PorousFlowBasicTHM<<<{"href": "../../../../syntax/PorousFlowBasicTHM/index.html"}>>>]
  dictator_name<<<{"description": "The name of the dictator user object that is created by this Action"}>>> = dictator
  add_darcy_aux<<<{"description": "Add AuxVariables that record Darcy velocity"}>>> = false
  fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = simple_fluid
  gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  multiply_by_density<<<{"description": "If true, then the Kernels for fluid flow are multiplied by the fluid density.  If false, this multiplication is not performed, which means the problem linearises, but that care must be taken when using other PorousFlow objects."}>>> = false
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = pp
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 0.001
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.05
  []
  [biot_mod]
    type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time.  Sometimes 1 / BiotModulus is called storativity", "href": "../../../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
    fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 2E9
    biot_coefficient<<<{"description": "Biot coefficient"}>>> = 1.0
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [sink]
    type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
    point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0 0 0'
    mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = -0.16E-3 # recall this is a volumetric flux because multiply_by_density = false in the Action, so this corresponds to a mass_flux of 0.16 kg/s/m because density=1000
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [pp]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 25
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
    end_point<<<{"description": "The ending of the line"}>>> = '100 0 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pp
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  dt = 200
  end_time = 1E3
  nl_abs_tol = 1e-10
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  [csv]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_on<<<{"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."}>>> = final
  []
[]

MOOSE agrees with the Theis solution:

Figure 2: Results of the fully-saturated, single-phase Theis simulation.

Flow from a source in a radial 1D problem admits a similarity solution, where the fluid pressure (and saturation, mass fraction, etc) is a function of the similarity variable $var element = document.getElementById("moose-equation-596ac53b-efb8-4ac8-9c17-2027decc5793");katex.render("\\zeta = r^2/t", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . This may be observed in Eq. (1). Using this fact, two-phase immiscible flow may be tested. The following input file describes injection of a gas phase into a fully liquid-saturated model at a constant rate:

# Two phase Theis problem: Flow from single source
# Constant rate injection 0.5 kg/s
# 1D cylindrical mesh

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 2000
  bias_x = 1.05
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [ppwater]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20e6
  []
  [sgas]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
  []
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [massfrac_ph0_sp0]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1
  []
  [massfrac_ph1_sp0]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
  []
  [flux0]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
  []
  [mass1]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = sgas
  []
  [flux1]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = sgas
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'ppwater sgas'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 2
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst<<<{"description": "Constant capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureConst.html"}>>>
    pc<<<{"description": "Constant capillary pressure (Pa). Default is 0"}>>> = 1e5
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid0]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
  [simple_fluid1]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 10
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-4
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow2PhasePS<<<{"description": "This Material calculates the 2 porepressures and the 2 saturations in a 2-phase situation, and derivatives of these with respect to the PorousFlowVariables.", "href": "../../../../source/materials/PorousFlow2PhasePS.html"}>>>
    phase0_porepressure<<<{"description": "Variable that is the porepressure of phase 0 (the liquid phase)"}>>> = ppwater
    phase1_saturation<<<{"description": "Variable that is the saturation of phase 1 (the gas phase)"}>>> = sgas
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
    mass_fraction_vars<<<{"description": "List of variables that represent the mass fractions.  Format is 'f_ph0^c0 f_ph0^c1 f_ph0^c2 ... f_ph0^c(N-2) f_ph1^c0 f_ph1^c1 fph1^c2 ... fph1^c(N-2) ... fphP^c0 f_phP^c1 fphP^c2 ... fphP^c(N-2)' where N=num_components and P=num_phases, and it is assumed that f_ph^c(N-1)=1-sum(f_ph^c,{c,0,N-2}) so that f_ph^c(N-1) need not be given.  If no variables are provided then num_phases=1=num_components."}>>> = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [simple_fluid0]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid0
    phase<<<{"description": "The phase number"}>>> = 0
    compute_enthalpy<<<{"description": "Compute the fluid enthalpy"}>>> = false
    compute_internal_energy<<<{"description": "Compute the fluid internal energy"}>>> = false
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid1
    phase<<<{"description": "The phase number"}>>> = 1
    compute_enthalpy<<<{"description": "Compute the fluid enthalpy"}>>> = false
    compute_internal_energy<<<{"description": "Compute the fluid internal energy"}>>> = false
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 1
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 1
    phase<<<{"description": "The phase number"}>>> = 1
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [rightwater]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "Value of the BC"}>>> = 20e6
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [source]
    type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
    point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0 0 0'
    mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = 0.5
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = sgas
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'gmres      asm      lu           NONZERO                   2               1E-8       1E-10 20'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  end_time = 1e4
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 10
    growth_factor = 2
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [line]
    type = NodalValueSampler<<<{"description": "Samples values of nodal variable(s).", "href": "../../../../source/vectorpostprocessors/NodalValueSampler.html"}>>>
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'ppwater sgas'
    execute_on<<<{"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."}>>> = 'timestep_end'
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [ppwater]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = ppwater
  []
  [sgas]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = sgas
  []
  [massgas]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component that this Postprocessor acts on"}>>> = 1
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = theis3
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  [csv]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_on<<<{"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."}>>> = timestep_end
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = final
  []
[]

Figure 3 shows the comparison of similarity solutions calculated with either fixed radial distance or fixed time. In this case, good agreement is observed between the two results for both liquid pressure and gas saturation.

Figure 3: Results of the 2-phase radial injection simulation.

Further tests of this kind, involving multi-component, multi-phase flow, including mutual dissolution of gas and liquid phases may be found here.

Peaceman borehole fluxes

The test suite that checks that the Peaceman flux $(2)var element = document.getElementById("moose-equation-6871a3c9-8e18-481d-b02e-fc687e4dfe33");katex.render("f(P_{i}, x_{i}) = W \\left|C(P_{i}-P_{\\mathrm{bh}})\\right| \\frac{k_{\\mathrm{r}}\\rho}{\\mu}(P_{i} - P_{\\mathrm{bh}}) ", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is correctly implemented. A vertical borehole is placed through the centre of a single element, and fluid flow to the borehole as a function of porepressure is measured. The tests are

A production borehole with $var element = document.getElementById("moose-equation-16708744-bca3-4016-8c3d-91f55513dde4");katex.render("P_{\\mathrm{bh}} = 0", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , with a fully-saturated medium.

# fully-saturated
# production
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1E7
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../../../source/userobjects/PorousFlowSumQuantity.html"}>>>
  []
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.5
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-7
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [bh]
    type = PorousFlowPeacemanBorehole<<<{"description": "Approximates a borehole in the mesh using the Peaceman approach, ie using a number of point sinks with given radii whose positions are read from a file.  NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../../../source/dirackernels/PorousFlowPeacemanBorehole.html"}>>>

    # 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<<<{"description": "The name of the variable that this residual object operates on"}>>> = 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<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = borehole_total_outflow_mass

    # The following file defines the polyline geometry
    # which is just two points in this particular example
    point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink.  The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii.  Each line in the file must contain a space-separated weight and coordinate, viz r x y z.  For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure.  If your file contains just one point, you must also specify the line_length and line_direction parameters.  Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = 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<<<{"description": "Modifying functions will be a function of either pressure and permeability (eg, for boreholes that pump fluids) or temperature and thermal conductivity (eg, for boreholes that pump pure heat with no fluid flow)"}>>> = pressure
    fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink.  For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0

    # The bottomhole pressure
    bottom_p_or_t<<<{"description": "For function_of=pressure, this function is the pressure at the bottom of the borehole, otherwise it is the temperature at the bottom of the borehole."}>>> = 0

    # In this example there is no increase of the wellbore pressure
    # due to gravity:
    unit_weight<<<{"description": "(fluid_density*gravitational_acceleration) as a vector pointing downwards.  Note that the borehole pressure at a given z position is bottom_p_or_t + unit_weight*(q - q_bottom), where q=(x,y,z) and q_bottom=(x,y,z) of the bottom point of the borehole.  The analogous formula holds for function_of=temperature.  If you don't want bottomhole pressure (or temperature) to vary in the borehole just set unit_weight=0.  Typical value is = (0,0,-1E4), for water"}>>> = '0 0 0'

    # PeacemanBoreholes should almost always have use_mobility = true
    use_mobility<<<{"description": "Multiply the flux by the fluid mobility"}>>> = true

    # This is a production wellbore (a sink of fluid that removes fluid from porespace)
    character<<<{"description": "If zero then borehole does nothing.  If positive the borehole acts as a sink (production well) for porepressure > borehole pressure, and does nothing otherwise.  If negative the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise.  The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1, but you can specify other quantities to provide an overall scaling to the flow if you like."}>>> = 1
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [bh_report]
    type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
    uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = borehole_total_outflow_mass
  []
  [fluid_mass0]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_begin
  []

  [fluid_mass1]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_end
  []

  [zmass_error]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = mass_bal_fcn
    execute_on<<<{"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."}>>> = timestep_end
    indirect_dependencies<<<{"description": "If the evaluated function depends on other postprocessors they must be listed here to ensure proper dependency resolution"}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []

  [p0]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [mass_bal_fcn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = abs((a-c+d)/2/(a+c))
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'a c d'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [usual]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = bh02
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"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."}>>> = timestep_end
[]

An injection borehole with $var element = document.getElementById("moose-equation-94831a45-6d8f-4b34-97a6-fb579c8de340");katex.render("P_{\\mathrm{bh}} = 10\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ MPa, with a fully-saturated medium.

# fully-saturated
# injection
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../../../source/userobjects/PorousFlowSumQuantity.html"}>>>
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.5
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-7
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [bh]
    type = PorousFlowPeacemanBorehole<<<{"description": "Approximates a borehole in the mesh using the Peaceman approach, ie using a number of point sinks with given radii whose positions are read from a file.  NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../../../source/dirackernels/PorousFlowPeacemanBorehole.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = borehole_total_outflow_mass
    point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink.  The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii.  Each line in the file must contain a space-separated weight and coordinate, viz r x y z.  For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure.  If your file contains just one point, you must also specify the line_length and line_direction parameters.  Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = bh03.bh
    function_of<<<{"description": "Modifying functions will be a function of either pressure and permeability (eg, for boreholes that pump fluids) or temperature and thermal conductivity (eg, for boreholes that pump pure heat with no fluid flow)"}>>> = pressure
    fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink.  For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0
    bottom_p_or_t<<<{"description": "For function_of=pressure, this function is the pressure at the bottom of the borehole, otherwise it is the temperature at the bottom of the borehole."}>>> = 'insitu_pp'
    unit_weight<<<{"description": "(fluid_density*gravitational_acceleration) as a vector pointing downwards.  Note that the borehole pressure at a given z position is bottom_p_or_t + unit_weight*(q - q_bottom), where q=(x,y,z) and q_bottom=(x,y,z) of the bottom point of the borehole.  The analogous formula holds for function_of=temperature.  If you don't want bottomhole pressure (or temperature) to vary in the borehole just set unit_weight=0.  Typical value is = (0,0,-1E4), for water"}>>> = '0 0 0'
    use_mobility<<<{"description": "Multiply the flux by the fluid mobility"}>>> = true
    character<<<{"description": "If zero then borehole does nothing.  If positive the borehole acts as a sink (production well) for porepressure > borehole pressure, and does nothing otherwise.  If negative the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise.  The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1, but you can specify other quantities to provide an overall scaling to the flow if you like."}>>> = -1
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [bh_report]
    type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
    uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = borehole_total_outflow_mass
  []

  [fluid_mass0]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_begin
  []

  [fluid_mass1]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_end
  []

  [zmass_error]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = mass_bal_fcn
    execute_on<<<{"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."}>>> = timestep_end
    indirect_dependencies<<<{"description": "If the evaluated function depends on other postprocessors they must be listed here to ensure proper dependency resolution"}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []

  [p0]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [mass_bal_fcn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = abs((a-c+d)/2/(a+c))
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'a c d'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []
  [insitu_pp]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '-2e7*z' #bh_bottom is located at z=-0.5
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [usual]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  end_time = 0.5
  dt = 1E-2
  solve_type = NEWTON
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = bh03
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"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."}>>> = timestep_end
[]

A production borehole with $var element = document.getElementById("moose-equation-8561b9b1-a32c-4478-b91a-9d9eab7b2151");katex.render("P_{\\mathrm{bh}} = -10\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ MPa with a fully-saturated medium, with use_mobility=true

# fully-saturated
# production
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [dts]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../../../source/functions/PiecewiseLinear.html"}>>>
    y<<<{"description": "The ordinate values"}>>> = '1E-2 1E-1 1 1E1 1E2 1E3'
    x<<<{"description": "The abscissa values"}>>> = '0 1E-1 1 1E1 1E2 1E3'
  []
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../../../source/userobjects/PorousFlowSumQuantity.html"}>>>
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC<<<{"description": "This Material calculates relative permeability of a phase using a model inspired by FLAC", "href": "../../../../source/materials/PorousFlowRelativePermeabilityFLAC.html"}>>>
    m<<<{"description": "relperm = (1 + m)seff^m - m seff^(m+1)"}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [bh]
    type = PorousFlowPeacemanBorehole<<<{"description": "Approximates a borehole in the mesh using the Peaceman approach, ie using a number of point sinks with given radii whose positions are read from a file.  NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../../../source/dirackernels/PorousFlowPeacemanBorehole.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = borehole_total_outflow_mass
    point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink.  The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii.  Each line in the file must contain a space-separated weight and coordinate, viz r x y z.  For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure.  If your file contains just one point, you must also specify the line_length and line_direction parameters.  Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = bh02.bh
    fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink.  For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0
    bottom_p_or_t<<<{"description": "For function_of=pressure, this function is the pressure at the bottom of the borehole, otherwise it is the temperature at the bottom of the borehole."}>>> = -1E6
    unit_weight<<<{"description": "(fluid_density*gravitational_acceleration) as a vector pointing downwards.  Note that the borehole pressure at a given z position is bottom_p_or_t + unit_weight*(q - q_bottom), where q=(x,y,z) and q_bottom=(x,y,z) of the bottom point of the borehole.  The analogous formula holds for function_of=temperature.  If you don't want bottomhole pressure (or temperature) to vary in the borehole just set unit_weight=0.  Typical value is = (0,0,-1E4), for water"}>>> = '0 0 0'
    use_mobility<<<{"description": "Multiply the flux by the fluid mobility"}>>> = true
    character<<<{"description": "If zero then borehole does nothing.  If positive the borehole acts as a sink (production well) for porepressure > borehole pressure, and does nothing otherwise.  If negative the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise.  The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1, but you can specify other quantities to provide an overall scaling to the flow if you like."}>>> = 1
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [bh_report]
    type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
    uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = borehole_total_outflow_mass
  []

  [fluid_mass0]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_begin
  []

  [fluid_mass1]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_end
  []

  [zmass_error]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = mass_bal_fcn
    execute_on<<<{"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."}>>> = timestep_end
    indirect_dependencies<<<{"description": "If the evaluated function depends on other postprocessors they must be listed here to ensure proper dependency resolution"}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []

  [p0]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
[]

[Functions]
  [mass_bal_fcn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = abs((a-c+d)/2/(a+c))
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'a c d'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [usual]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  end_time = 1E3
  solve_type = NEWTON
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = dts
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = bh04
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"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."}>>> = timestep_end
[]

An injection borehole with $var element = document.getElementById("moose-equation-2249c6cb-c02b-4ff1-bea8-2c959b03fd76");katex.render("P_{\\mathrm{bh}} = 0", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ with an unsaturated medium, with use_mobility=true

# unsaturated
# injection
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  zmin = -1
  zmax = 1
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [dts]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../../../source/functions/PiecewiseLinear.html"}>>>
    y<<<{"description": "The ordinate values"}>>> = '500 500 1E1'
    x<<<{"description": "The abscissa values"}>>> = '4000 5000 6500'
  []
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = -2E5
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../../../source/userobjects/PorousFlowSumQuantity.html"}>>>
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC<<<{"description": "This Material calculates relative permeability of a phase using a model inspired by FLAC", "href": "../../../../source/materials/PorousFlowRelativePermeabilityFLAC.html"}>>>
    m<<<{"description": "relperm = (1 + m)seff^m - m seff^(m+1)"}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [bh]
    type = PorousFlowPeacemanBorehole<<<{"description": "Approximates a borehole in the mesh using the Peaceman approach, ie using a number of point sinks with given radii whose positions are read from a file.  NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../../../source/dirackernels/PorousFlowPeacemanBorehole.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = borehole_total_outflow_mass
    point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink.  The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii.  Each line in the file must contain a space-separated weight and coordinate, viz r x y z.  For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure.  If your file contains just one point, you must also specify the line_length and line_direction parameters.  Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = bh03.bh
    fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink.  For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0
    bottom_p_or_t<<<{"description": "For function_of=pressure, this function is the pressure at the bottom of the borehole, otherwise it is the temperature at the bottom of the borehole."}>>> = 0
    unit_weight<<<{"description": "(fluid_density*gravitational_acceleration) as a vector pointing downwards.  Note that the borehole pressure at a given z position is bottom_p_or_t + unit_weight*(q - q_bottom), where q=(x,y,z) and q_bottom=(x,y,z) of the bottom point of the borehole.  The analogous formula holds for function_of=temperature.  If you don't want bottomhole pressure (or temperature) to vary in the borehole just set unit_weight=0.  Typical value is = (0,0,-1E4), for water"}>>> = '0 0 0'
    use_mobility<<<{"description": "Multiply the flux by the fluid mobility"}>>> = true
    character<<<{"description": "If zero then borehole does nothing.  If positive the borehole acts as a sink (production well) for porepressure > borehole pressure, and does nothing otherwise.  If negative the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise.  The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1, but you can specify other quantities to provide an overall scaling to the flow if you like."}>>> = -1
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [bh_report]
    type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
    uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = borehole_total_outflow_mass
  []

  [fluid_mass0]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_begin
  []

  [fluid_mass1]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = timestep_end
  []

  [zmass_error]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = mass_bal_fcn
    execute_on<<<{"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."}>>> = timestep_end
    indirect_dependencies<<<{"description": "If the evaluated function depends on other postprocessors they must be listed here to ensure proper dependency resolution"}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []

  [p0]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pp
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    execute_on<<<{"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."}>>> = timestep_end
  []
[]

[Functions]
  [mass_bal_fcn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = abs((a-c+d)/2/(a+c))
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'a c d'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'fluid_mass1 fluid_mass0 bh_report'
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [usual]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  end_time = 6500
  solve_type = NEWTON
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = dts
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = bh05
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = false
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"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."}>>> = timestep_end
[]

Further commentary, and the results demonstrating that MOOSE is correct may be found in the page discussing the theory behind Dirac Kernels in PorousFlow.

Comparison with a steady-state 2D analytic solution

The test

# Comparison with analytical solution for cylindrically-symmetric situation
[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = FileMesh
  file = bh07_input.e
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [dts]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../../../source/functions/PiecewiseLinear.html"}>>>
    y<<<{"description": "The ordinate values"}>>> = '1000 10000'
    x<<<{"description": "The abscissa values"}>>> = '100 1000'
  []
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1E7
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
  []
  [fflux]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [fix_outer]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = perimeter
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    value<<<{"description": "Value of the BC"}>>> = 1E7
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pp'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 1
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 1
  []
  [borehole_total_outflow_mass]
    type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../../../source/userobjects/PorousFlowSumQuantity.html"}>>>
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.8
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2e9
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
  []
  [ppss]
    type = PorousFlow1PhaseP<<<{"description": "This Material is used for the partially saturated single-phase situation where porepressure is the primary variable", "href": "../../../../source/materials/PorousFlow1PhaseP.html"}>>>
    porepressure<<<{"description": "Variable that represents the porepressure of the single phase"}>>> = pp
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
  [massfrac]
    type = PorousFlowMassFraction<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions out of the individual mass fractions", "href": "../../../../source/materials/PorousFlowMassFraction.html"}>>>
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid<<<{"description": "This Material calculates fluid properties at the quadpoints or nodes for a single component fluid", "href": "../../../../source/materials/PorousFlowSingleComponentFluid.html"}>>>
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
    phase<<<{"description": "The phase number"}>>> = 0
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-11 0 0 0 1E-11 0 0 0 1E-11'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityFLAC<<<{"description": "This Material calculates relative permeability of a phase using a model inspired by FLAC", "href": "../../../../source/materials/PorousFlowRelativePermeabilityFLAC.html"}>>>
    m<<<{"description": "relperm = (1 + m)seff^m - m seff^(m+1)"}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [bh]
    type = PorousFlowPeacemanBorehole<<<{"description": "Approximates a borehole in the mesh using the Peaceman approach, ie using a number of point sinks with given radii whose positions are read from a file.  NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../../../source/dirackernels/PorousFlowPeacemanBorehole.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = borehole_total_outflow_mass
    point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink.  The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii.  Each line in the file must contain a space-separated weight and coordinate, viz r x y z.  For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure.  If your file contains just one point, you must also specify the line_length and line_direction parameters.  Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = bh07.bh
    fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink.  For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0
    bottom_p_or_t<<<{"description": "For function_of=pressure, this function is the pressure at the bottom of the borehole, otherwise it is the temperature at the bottom of the borehole."}>>> = 0
    unit_weight<<<{"description": "(fluid_density*gravitational_acceleration) as a vector pointing downwards.  Note that the borehole pressure at a given z position is bottom_p_or_t + unit_weight*(q - q_bottom), where q=(x,y,z) and q_bottom=(x,y,z) of the bottom point of the borehole.  The analogous formula holds for function_of=temperature.  If you don't want bottomhole pressure (or temperature) to vary in the borehole just set unit_weight=0.  Typical value is = (0,0,-1E4), for water"}>>> = '0 0 0'
    use_mobility<<<{"description": "Multiply the flux by the fluid mobility"}>>> = true
    re_constant<<<{"description": "The dimensionless constant used in evaluating the borehole effective radius.  This depends on the meshing scheme.  Peacemann finite-difference calculations give 0.28, while for rectangular finite elements the result is closer to 0.1594.  (See  Eqn(4.13) of Z Chen, Y Zhang, Well flow models for various numerical methods, Int J Num Analysis and Modeling, 3 (2008) 375-388.)"}>>> = 0.1594  # use Chen and Zhang version
    character<<<{"description": "If zero then borehole does nothing.  If positive the borehole acts as a sink (production well) for porepressure > borehole pressure, and does nothing otherwise.  If negative the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise.  The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1, but you can specify other quantities to provide an overall scaling to the flow if you like."}>>> = 2 # double the strength because bh07.bh only fills half the mesh
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [bh_report]
    type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
    uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = borehole_total_outflow_mass
    execute_on<<<{"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."}>>> = 'initial timestep_end'
  []
  [fluid_mass]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    execute_on<<<{"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."}>>> = 'initial timestep_end'
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [pp]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pp
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
    end_point<<<{"description": "The ending of the line"}>>> = '300 0 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 300
    execute_on<<<{"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."}>>> = timestep_end
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [usual]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-10 1E-10 10000 30'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  end_time = 1E3
  solve_type = NEWTON
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    # get only marginally better results for smaller time steps
    type = FunctionDT
    function = dts
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = bh07
  [along_line]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    execute_on<<<{"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."}>>> = final
  []
  [exodus]
    type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../../../../source/outputs/Exodus.html"}>>>
    execute_on<<<{"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."}>>> = 'initial final'
  []
[]

checks whether the Peaceman borehole extracts fluid at the correct rate and the corresponding reduction in porepressure agrees with an analytic solution of Laplace's equation. Details and results demonstrating the correctness of PorousFlow can be found in the page discussing the theory behind Dirac Kernels in PorousFlow.

(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/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/squarepulse1.i)

# Test PorousFlowSquarePulsePointSource DiracKernel

[Mesh]
  type = GeneratedMesh
  dim = 2
  bias_x = 1.1
  bias_y = 1.1
  ymax = 1
  xmax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
  []
[]

[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
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    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.2
  []
[]

[Postprocessors]
  [total_mass]
    type = PorousFlowFluidMass
    execute_on = 'initial timestep_end'
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  nl_abs_tol = 1e-14
  dt = 200
  end_time = 2000
[]

[Outputs]
  perf_graph = true
  file_base = squarepulse1
  csv = true
  execute_on = 'initial timestep_end'
  [con]
    output_linear = true
    type = Console
  []
[]

[ICs]
  [PressureIC]
    variable = pp
    type = ConstantIC
    value = 20e6
  []
[]

[DiracKernels]
  [sink1]
    type = PorousFlowSquarePulsePointSource
    start_time = 100
    end_time = 300
    point = '0.5 0.5 0'
    mass_flux = -0.1
    variable = pp
  []
  [sink]
    type = PorousFlowSquarePulsePointSource
    start_time = 600
    end_time = 1400
    point = '0.5 0.5 0'
    mass_flux = -0.1
    variable = pp
  []
  [source]
    point = '0.5 0.5 0'
    start_time = 1500
    mass_flux = 0.2
    end_time = 2000
    variable = pp
    type = PorousFlowSquarePulsePointSource
  []
[]

(modules/porous_flow/test/tests/dirackernels/theis_rz.i)

# Theis problem: Flow to single sink using BasicTHM
# SinglePhase
# RZ mesh

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 20
  xmax = 100
  bias_x = 1.05
  coord_type = RZ
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    initial_condition = 20E6
  []
[]

[PorousFlowBasicTHM]
  dictator_name = dictator
  add_darcy_aux = false
  fp = simple_fluid
  gravity = '0 0 0'
  multiply_by_density = false
  porepressure = pp
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    viscosity = 0.001
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.05
  []
  [biot_mod]
    type = PorousFlowConstantBiotModulus
    fluid_bulk_modulus = 2E9
    biot_coefficient = 1.0
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
  []
[]

[DiracKernels]
  [sink]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = -0.16E-3 # recall this is a volumetric flux because multiply_by_density = false in the Action, so this corresponds to a mass_flux of 0.16 kg/s/m because density=1000
    variable = pp
  []
[]

[VectorPostprocessors]
  [pp]
    type = LineValueSampler
    num_points = 25
    start_point = '0 0 0'
    end_point = '100 0 0'
    sort_by = x
    variable = pp
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 200
  end_time = 1E3
  nl_abs_tol = 1e-10
[]

[Outputs]
  perf_graph = true
  [csv]
    type = CSV
    execute_on = final
  []
[]

(modules/porous_flow/test/tests/dirackernels/theis3.i)

# Two phase Theis problem: Flow from single source
# Constant rate injection 0.5 kg/s
# 1D cylindrical mesh

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 2000
  bias_x = 1.05
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [ppwater]
    initial_condition = 20e6
  []
  [sgas]
    initial_condition = 0
  []
[]

[AuxVariables]
  [massfrac_ph0_sp0]
    initial_condition = 1
  []
  [massfrac_ph1_sp0]
    initial_condition = 0
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = ppwater
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = ppwater
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = sgas
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = sgas
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'ppwater sgas'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 1e5
  []
[]

[FluidProperties]
  [simple_fluid0]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 1000
    viscosity = 1e-3
    thermal_expansion = 0
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 10
    viscosity = 1e-4
    thermal_expansion = 0
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow2PhasePS
    phase0_porepressure = ppwater
    phase1_saturation = sgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [simple_fluid0]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
    compute_enthalpy = false
    compute_internal_energy = false
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
    compute_enthalpy = false
    compute_internal_energy = false
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 1
  []
[]

[BCs]
  [rightwater]
    type = DirichletBC
    boundary = right
    value = 20e6
    variable = ppwater
  []
[]

[DiracKernels]
  [source]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 0.5
    variable = sgas
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-8       1E-10 20'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1e4
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 10
    growth_factor = 2
  []
[]

[VectorPostprocessors]
  [line]
    type = NodalValueSampler
    sort_by = x
    variable = 'ppwater sgas'
    execute_on = 'timestep_end'
  []
[]

[Postprocessors]
  [ppwater]
    type = PointValue
    point = '4 0 0'
    variable = ppwater
  []
  [sgas]
    type = PointValue
    point = '4 0 0'
    variable = sgas
  []
  [massgas]
    type = PorousFlowFluidMass
    fluid_component = 1
  []
[]

[Outputs]
  file_base = theis3
  print_linear_residuals = false
  perf_graph = true
  [csv]
    type = CSV
    execute_on = timestep_end
    execute_vector_postprocessors_on = final
  []
[]

(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/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/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/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/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'
  []
[]

Geometric tests
Basic point sources / sinks
Theis tests
Peaceman borehole fluxes
Comparison with a steady - state 2 D analytic solution