FVPorousFlowAdvectiveFlux

Advective Darcy flux

This FVKernel implements the strong form of $var element = document.getElementById("moose-equation-cb6e1d03-7874-4080-9102-424dcdc819fe");katex.render(" -\\nabla\\cdot \\sum_{\\beta}\\chi_{\\beta}^{\\kappa} \\rho_{\\beta}\\frac{k\\,k_{\\mathrm{r,}\\beta}}{\\mu_{\\beta}}(\\nabla P_{\\beta} - \\rho_{\\beta} \\mathbf{g})", 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 all parameters are defined in the nomenclature.

Input Parameters

PorousFlowDictatorThe PorousFlowDictator UserObject
C++ Type:UserObjectName
Controllable:No
Description:The PorousFlowDictator UserObject
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
fluid_component0The fluid component for this kernel
Default:0
C++ Type:unsigned int
Controllable:No
Description:The fluid component for this kernel
gravity0 0 -9.81Gravity vector. Defaults to (0, 0, -9.81)
Default:0 0 -9.81
C++ Type:libMesh::VectorValue<double>
Unit:(no unit assumed)
Controllable:No
Description:Gravity vector. Defaults to (0, 0, -9.81)
matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)

Optional Parameters

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

Contribution To Tagged Field Data Parameters

boundaries_to_avoidThe set of sidesets to not execute this FVFluxKernel on. This takes precedence over force_boundary_execution to restrict to less external boundaries. By default flux kernels are executed on all internal boundaries and Dirichlet boundary conditions.
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The set of sidesets to not execute this FVFluxKernel on. This takes precedence over force_boundary_execution to restrict to less external boundaries. By default flux kernels are executed on all internal boundaries and Dirichlet boundary conditions.
boundaries_to_forceThe set of sidesets to force execution of this FVFluxKernel on. Setting force_boundary_execution to true is equivalent to listing all external mesh boundaries in this parameter.
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The set of sidesets to force execution of this FVFluxKernel on. Setting force_boundary_execution to true is equivalent to listing all external mesh boundaries in this parameter.
force_boundary_executionFalseWhether to force execution of this object on all external boundaries.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to force execution of this object on all external boundaries.

Boundary Execution Modification Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

ghost_layers2The number of layers of elements to ghost.
Default:2
C++ Type:unsigned short
Controllable:No
Description:The number of layers of elements to ghost.
use_point_neighborsFalseWhether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.

Parallel Ghosting Parameters

prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

Input Files

(modules/porous_flow/test/tests/poroperm/PermTensorFromVar01_fv.i)
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_2phase_fv.i)
(modules/porous_flow/test/tests/gravity/grav02b_fv.i)
(modules/porous_flow/test/tests/jacobian/fv_mass_flux.i)
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_fv.i)
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_fully_saturated_fv.i)
(modules/porous_flow/test/tests/poroperm/PermFromPoro01_fv.i)
(modules/porous_flow/test/tests/heat_advection/heat_advection_1d_fv.i)
(modules/porous_flow/test/tests/heterogeneous_materials/constant_poroperm_fv.i)
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_2phasePS_fv.i)
(modules/porous_flow/test/tests/gravity/grav02e_fv.i)
(modules/porous_flow/test/tests/gravity/grav01a_fv.i)
(modules/porous_flow/test/tests/dispersion/disp01_fv.i)
(modules/porous_flow/test/tests/poroperm/PermFromPoro03_fv.i)
(modules/porous_flow/examples/fluidflower/fluidflower.i)

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

# Testing permeability calculated from scalar and tensor
# Trivial test, checking calculated permeability is correct
# k = k_anisotropy * perm

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 3
    xmin = 0
    xmax = 3
  []
[]

[GlobalParams]
  block = 0
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    type = MooseVariableFVReal
    [FVInitialCondition]
      type = FVConstantIC
      value = 0
    []
  []
[]

[FVKernels]
  [flux]
    type = FVPorousFlowAdvectiveFlux
    gravity = '0 0 0'
    variable = pp
  []
[]

[FVBCs]
  [ptop]
    type = FVDirichletBC
    variable = pp
    boundary = right
    value = 0
  []
  [pbase]
    type = FVDirichletBC
    variable = pp
    boundary = left
    value = 1
  []
[]

[AuxVariables]
  [perm_var]
    type = MooseVariableFVReal
  []
  [perm_x]
    type = MooseVariableFVReal
  []
  [perm_y]
    type = MooseVariableFVReal
  []
  [perm_z]
    type = MooseVariableFVReal
  []
[]

[AuxKernels]
  [perm_var]
    type = ConstantAux
    value = 2
    variable = perm_var
  []
  [perm_x]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_x
    row = 0
    column = 0
  []
  [perm_y]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_y
    row = 1
    column = 1
  []
  [perm_z]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_z
    row = 2
    column = 2
  []
[]

[Postprocessors]
  [perm_x_left]
    type = PointValue
    variable = perm_x
    point = '0.5 0 0'
  []
  [perm_y_left]
    type = PointValue
    variable = perm_y
    point = '0.5 0 0'
  []
  [perm_z_left]
    type = PointValue
    variable = perm_z
    point = '0.5 0 0'
  []
  [perm_x_right]
    type = PointValue
    variable = perm_x
    point = '2.5 0 0'
  []
  [perm_y_right]
    type = PointValue
    variable = perm_y
    point = '2.5 0 0'
  []
  [perm_z_right]
    type = PointValue
    variable = perm_z
    point = '2.5 0 0'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    # unimportant in this fully-saturated test
    m = 0.8
    alpha = 1e-4
  []
[]

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

[Materials]
  [permeability]
    type = ADPorousFlowPermeabilityTensorFromVar
    k_anisotropy = '1 0 0  0 2 0  0 0 0.1'
    perm = perm_var
  []
  [temperature]
    type = ADPorousFlowTemperature
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [eff_fluid_pressure]
    type = ADPorousFlowEffectiveFluidPressure
  []
  [ppss]
    type = ADPorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 0 # unimportant in this fully-saturated situation
    phase = 0
  []
[]

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

[Executioner]
  type = Steady
  solve_type = Newton
  l_tol = 1E-5
  nl_abs_tol = 1E-3
  nl_rel_tol = 1E-8
  l_max_its = 200
  nl_max_its = 400
[]

[Outputs]
  file_base = 'PermTensorFromVar01_out'
  csv = true
  execute_on = 'timestep_end'
[]

(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_2phase_fv.i)

# Pressure pulse in 1D with 2 phases (with one having zero saturation), 2components - transient using FV
[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 10
    xmin = 0
    xmax = 100
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [ppwater]
    type = MooseVariableFVReal
    initial_condition = 2E6
  []
  [ppgas]
    type = MooseVariableFVReal
    initial_condition = 2E6
  []
[]

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

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = ppwater
  []
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    variable = ppwater
    gravity = '0 0 0'
    fluid_component = 0
  []
  [mass1]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 1
    variable = ppgas
  []
  [flux1]
    type = FVPorousFlowAdvectiveFlux
    variable = ppgas
    gravity = '0 0 0'
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'ppwater ppgas'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
[]

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

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
  []
  [ppss]
    type = ADPorousFlow2PhasePP
    phase0_porepressure = ppwater
    phase1_porepressure = ppgas
    capillary_pressure = pc
  []
  [massfrac]
    type = ADPorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [simple_fluid0]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm_water]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
  [relperm_gas]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 1
    phase = 1
  []
[]

[FVBCs]
  [leftwater]
    type = FVDirichletBC
    boundary = left
    value = 3E6
    variable = ppwater
  []
  [leftgas]
    type = FVDirichletBC
    boundary = left
    value = 3E6
    variable = ppgas
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2               1E-12'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E3
  end_time = 1E4
[]

[Postprocessors]
  [p005]
    type = PointValue
    variable = ppwater
    point = '5 0 0'
    execute_on = 'initial timestep_end'
  []
  [p015]
    type = PointValue
    variable = ppwater
    point = '15 0 0'
    execute_on = 'initial timestep_end'
  []
  [p025]
    type = PointValue
    variable = ppwater
    point = '25 0 0'
    execute_on = 'initial timestep_end'
  []
  [p035]
    type = PointValue
    variable = ppwater
    point = '35 0 0'
    execute_on = 'initial timestep_end'
  []
  [p045]
    type = PointValue
    variable = ppwater
    point = '45 0 0'
    execute_on = 'initial timestep_end'
  []
  [p055]
    type = PointValue
    variable = ppwater
    point = '55 0 0'
    execute_on = 'initial timestep_end'
  []
  [p065]
    type = PointValue
    variable = ppwater
    point = '65 0 0'
    execute_on = 'initial timestep_end'
  []
  [p075]
    type = PointValue
    variable = ppwater
    point = '75 0 0'
    execute_on = 'initial timestep_end'
  []
  [p085]
    type = PointValue
    variable = ppwater
    point = '85 0 0'
    execute_on = 'initial timestep_end'
  []
  [p095]
    type = PointValue
    variable = ppwater
    point = '95 0 0'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = pressure_pulse_1d_2phase_fv
  print_linear_residuals = false
  csv = true
[]

(modules/porous_flow/test/tests/gravity/grav02b_fv.i)

# Checking that gravity head is established in the steady-state situation when 0<saturation<1 (note the strictly less-than).
# 2phase (PP), 2components, vanGenuchten, constant fluid bulk-moduli for each phase, constant viscosity, constant permeability, Corey relative perm
# For better agreement with the analytical solution (ana_pp), just increase nx

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmin = -1
  xmax = 0
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [ppwater]
    type = MooseVariableFVReal
    initial_condition = -1.0
  []
  [ppgas]
    type = MooseVariableFVReal
    initial_condition = 0
  []
[]

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

[FVKernels]
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 0
    variable = ppwater
    gravity = '-1 0 0'
  []
  [flux1]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 1
    variable = ppgas
    gravity = '-1 0 0'
  []
[]

[FVBCs]
  [ppwater]
    type = FVDirichletBC
    boundary = right
    variable = ppwater
    value = -1
  []
  [ppgas]
    type = FVDirichletBC
    boundary = right
    variable = ppgas
    value = 0
  []
[]

[Functions]
  [ana_ppwater]
    type = ParsedFunction
    symbol_names = 'g B p0 rho0'
    symbol_values = '1 2 pp_water_top 1'
    expression = '-B*log(exp(-p0/B)+g*rho0*x/B)' # expected pp at base
  []
  [ana_ppgas]
    type = ParsedFunction
    symbol_names = 'g B p0 rho0'
    symbol_values = '1 1 pp_gas_top 0.1'
    expression = '-B*log(exp(-p0/B)+g*rho0*x/B)' # expected pp at base
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'ppwater ppgas'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  []
[]

[FluidProperties]
  [simple_fluid0]
    type = SimpleFluidProperties
    bulk_modulus = 2
    density0 = 1
    viscosity = 1
    thermal_expansion = 0
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 1
    density0 = 0.1
    viscosity = 0.5
    thermal_expansion = 0
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
  []
  [ppss]
    type = ADPorousFlow2PhasePP
    phase0_porepressure = ppwater
    phase1_porepressure = ppgas
    capillary_pressure = pc
  []
  [massfrac]
    type = ADPorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [simple_fluid0]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1 0 0  0 2 0  0 0 3'
  []
  [relperm_water]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
  [relperm_gas]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 1
    phase = 1
  []
[]

[Postprocessors]
  [pp_water_top]
    type = PointValue
    variable = ppwater
    point = '0 0 0'
  []
  [pp_water_base]
    type = PointValue
    variable = ppwater
    point = '-1 0 0'
  []
  [pp_water_analytical]
    type = FunctionValuePostprocessor
    function = ana_ppwater
    point = '-1 0 0'
  []
  [pp_gas_top]
    type = PointValue
    variable = ppgas
    point = '0 0 0'
  []
  [pp_gas_base]
    type = PointValue
    variable = ppgas
    point = '-1 0 0'
  []
  [pp_gas_analytical]
    type = FunctionValuePostprocessor
    function = ana_ppgas
    point = '-1 0 0'
  []
[]

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

[Executioner]
  type = Steady
  solve_type = Newton
[]

[Outputs]
  [csv]
    type = CSV
  []
[]

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

# Verify Jacobian of FV advective flux and mass time derivative kernels

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

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    family = MONOMIAL
    order = CONSTANT
    fv = true
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [flux]
    type = FVPorousFlowAdvectiveFlux
    variable = pp
    gravity = '0 -10 0'
    fluid_component = 0
  []
[]

[ICs]
  [pp]
    type = RandomIC
    variable = pp
    min = 1e6
    max = 2e6
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-5
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 1000
    thermal_expansion = 0
    viscosity = 1e-3
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
    temperature = 293
  []
  [ppss]
    type = ADPorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityConst
    kr = 1
    phase = 0
  []
[]

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

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

(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_fv.i)

# Pressure pulse in 1D with 1 phase - transient FV model

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmin = 0
  xmax = 100
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 2E6
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [flux]
    type = FVPorousFlowAdvectiveFlux
    variable = pp
    gravity = '0 0 0'
    fluid_component = 0
  []
[]

[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
    viscosity = 1e-3
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
    temperature = 293
  []
  [ppss]
    type = ADPorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityConst
    kr = 1
    phase = 0
  []
[]

[FVBCs]
  [left]
    type = FVPorousFlowAdvectiveFluxBC
    boundary = left
    porepressure_value = 3E6
    variable = pp
    gravity = '0 0 0'
    fluid_component = 0
    phase = 0
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E3
  end_time = 1E4
[]

[Postprocessors]
  [p005]
    type = PointValue
    variable = pp
    point = '5 0 0'
    execute_on = 'initial timestep_end'
  []
  [p015]
    type = PointValue
    variable = pp
    point = '15 0 0'
    execute_on = 'initial timestep_end'
  []
  [p025]
    type = PointValue
    variable = pp
    point = '25 0 0'
    execute_on = 'initial timestep_end'
  []
  [p035]
    type = PointValue
    variable = pp
    point = '35 0 0'
    execute_on = 'initial timestep_end'
  []
  [p045]
    type = PointValue
    variable = pp
    point = '45 0 0'
    execute_on = 'initial timestep_end'
  []
  [p055]
    type = PointValue
    variable = pp
    point = '55 0 0'
    execute_on = 'initial timestep_end'
  []
  [p065]
    type = PointValue
    variable = pp
    point = '65 0 0'
    execute_on = 'initial timestep_end'
  []
  [p075]
    type = PointValue
    variable = pp
    point = '75 0 0'
    execute_on = 'initial timestep_end'
  []
  [p085]
    type = PointValue
    variable = pp
    point = '85 0 0'
    execute_on = 'initial timestep_end'
  []
  [p095]
    type = PointValue
    variable = pp
    point = '95 0 0'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = pressure_pulse_1d_fv
  print_linear_residuals = false
  csv = true
[]

(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_fully_saturated_fv.i)

# Pressure pulse in 1D with 1 phase fully saturated - transient FV model

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmin = 0
  xmax = 100
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    type = MooseVariableFVReal
    initial_condition = 2E6
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [flux]
    type = FVPorousFlowAdvectiveFlux
    variable = pp
    gravity = '0 0 0'
    fluid_component = 0
  []
[]

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

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 1000
    thermal_expansion = 0
    viscosity = 1e-3
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
    temperature = 293
  []
  [ppss]
    type = ADPorousFlow1PhaseFullySaturated
    porepressure = pp
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityConst
    kr = 1
    phase = 0
  []
[]

[FVBCs]
  [left]
    type = FVPorousFlowAdvectiveFluxBC
    boundary = left
    porepressure_value = 3E6
    variable = pp
    gravity = '0 0 0'
    fluid_component = 0
    phase = 0
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E3
  end_time = 1E4
[]

[Postprocessors]
  [p005]
    type = PointValue
    variable = pp
    point = '5 0 0'
    execute_on = 'initial timestep_end'
  []
  [p015]
    type = PointValue
    variable = pp
    point = '15 0 0'
    execute_on = 'initial timestep_end'
  []
  [p025]
    type = PointValue
    variable = pp
    point = '25 0 0'
    execute_on = 'initial timestep_end'
  []
  [p035]
    type = PointValue
    variable = pp
    point = '35 0 0'
    execute_on = 'initial timestep_end'
  []
  [p045]
    type = PointValue
    variable = pp
    point = '45 0 0'
    execute_on = 'initial timestep_end'
  []
  [p055]
    type = PointValue
    variable = pp
    point = '55 0 0'
    execute_on = 'initial timestep_end'
  []
  [p065]
    type = PointValue
    variable = pp
    point = '65 0 0'
    execute_on = 'initial timestep_end'
  []
  [p075]
    type = PointValue
    variable = pp
    point = '75 0 0'
    execute_on = 'initial timestep_end'
  []
  [p085]
    type = PointValue
    variable = pp
    point = '85 0 0'
    execute_on = 'initial timestep_end'
  []
  [p095]
    type = PointValue
    variable = pp
    point = '95 0 0'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = pressure_pulse_1d_fv
  print_linear_residuals = false
  csv = true
[]

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

# Testing permeability from porosity
# Trivial test, checking calculated permeability is correct
# k = k_anisotropic * f * d^2 * phi^n / (1-phi)^m

[Mesh]
  [mesh]
  type = GeneratedMeshGenerator
    dim = 1
    nx = 3
    xmin = 0
    xmax = 3
  []
[]

[GlobalParams]
  block = 0
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    type = MooseVariableFVReal
    [FVInitialCondition]
      type = FVConstantIC
      value = 0
    []
  []
[]

[FVKernels]
  [flux]
    type = FVPorousFlowAdvectiveFlux
    gravity = '0 0 0'
    variable = pp
  []
[]

[FVBCs]
  [ptop]
    type = FVDirichletBC
    variable = pp
    boundary = right
    value = 0
  []
  [pbase]
    type = FVDirichletBC
    variable = pp
    boundary = left
    value = 1
  []
[]

[AuxVariables]
  [poro]
    type = MooseVariableFVReal
  []
  [perm_x]
    type = MooseVariableFVReal
  []
  [perm_y]
    type = MooseVariableFVReal
  []
  [perm_z]
    type = MooseVariableFVReal
  []
[]

[AuxKernels]
  [poro]
    type = ADPorousFlowPropertyAux
    property = porosity
    variable = poro
  []
  [perm_x]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_x
    row = 0
    column = 0
  []
  [perm_y]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_y
    row = 1
    column = 1
  []
  [perm_z]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_z
    row = 2
    column = 2
  []
[]

[Postprocessors]
  [perm_x_bottom]
    type = PointValue
    variable = perm_x
    point = '0 0 0'
  []
  [perm_y_bottom]
    type = PointValue
    variable = perm_y
    point = '0 0 0'
  []
  [perm_z_bottom]
    type = PointValue
    variable = perm_z
    point = '0 0 0'
  []
  [perm_x_top]
    type = PointValue
    variable = perm_x
    point = '3 0 0'
  []
  [perm_y_top]
    type = PointValue
    variable = perm_y
    point = '3 0 0'
  []
  [perm_z_top]
    type = PointValue
    variable = perm_z
    point = '3 0 0'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    # unimportant in this fully-saturated test
    m = 0.8
    alpha = 1e-4
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2.2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]

[Materials]
  [permeability]
    type = ADPorousFlowPermeabilityKozenyCarman
    k_anisotropy = '1 0 0  0 2 0  0 0 0.1'
    poroperm_function = kozeny_carman_fd2
    f = 0.1
    d = 5
    m = 2
    n = 7
  []
  [temperature]
    type = ADPorousFlowTemperature
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [eff_fluid_pressure]
    type = ADPorousFlowEffectiveFluidPressure
  []
  [ppss]
    type = ADPorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 0 # unimportant in this fully-saturated situation
    phase = 0
  []
[]

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

[Executioner]
  solve_type = Newton
  type = Steady
  l_tol = 1E-5
  nl_abs_tol = 1E-3
  nl_rel_tol = 1E-8
  l_max_its = 200
  nl_max_its = 400
[]

[Outputs]
  file_base = 'PermFromPoro01_out'
  csv = true
  execute_on = 'timestep_end'
[]

(modules/porous_flow/test/tests/heat_advection/heat_advection_1d_fv.i)

# 1phase, heat advecting with a moving fluid using FV

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 50
    xmin = 0
    xmax = 1
  []
[]

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

[Variables]
  [temp]
    type = MooseVariableFVReal
  []
  [pp]
    type = MooseVariableFVReal
  []
[]

[FVICs]
  [pp]
    type = FVFunctionIC
    variable = pp
    function = '1-x'
  []
  [temp]
    type = FVFunctionIC
    variable = temp
    function = 'if(x<0.02, 300, 200)'
  []
[]

[FVBCs]
  [pp0]
    type = FVDirichletBC
    variable = pp
    boundary = left
    value = 1
  []
  [pp1]
    type = FVDirichletBC
    variable = pp
    boundary = right
    value = 0
  []
  [hot]
    type = FVDirichletBC
    variable = temp
    boundary = left
    value = 300
  []
  [cold]
    type = FVDirichletBC
    variable = temp
    boundary = right
    value = 200
  []
[]

[FVKernels]
  [mass_dot]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [advection]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pp
  []
  [energy_dot]
    type = FVPorousFlowEnergyTimeDerivative
    variable = temp
  []
  [heat_advection]
    type = FVPorousFlowHeatAdvection
    variable = temp
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.6
    alpha = 1.3
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 100
    density0 = 1000
    viscosity = 4.4
    thermal_expansion = 0
    cv = 2
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
    temperature = temp
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.2
  []
  [rock_heat]
    type = ADPorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1.0
    density = 125
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 2 0 0 0 3'
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [PS]
    type = ADPorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 0.01
  end_time = 0.6
[]

[VectorPostprocessors]
  [T]
    type = ElementValueSampler
    sort_by = x
    variable = 'temp'
  []
[]

[Outputs]
  [csv]
    type = CSV
    execute_vector_postprocessors_on = final
  []
[]

(modules/porous_flow/test/tests/heterogeneous_materials/constant_poroperm_fv.i)

# Assign porosity and permeability variables from constant AuxVariables to create
# a heterogeneous model and solve with FV variables

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 3
    ny = 3
    nz = 3
    xmax = 3
    ymax = 3
    zmax = 3
  []
[]

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

[Variables]
  [ppwater]
    type = MooseVariableFVReal
    initial_condition = 1.5e6
  []
[]

[AuxVariables]
  [poro]
    type = MooseVariableFVReal
  []
  [permxx]
    type = MooseVariableFVReal
  []
  [permxy]
    type = MooseVariableFVReal
  []
  [permxz]
    type = MooseVariableFVReal
  []
  [permyx]
    type = MooseVariableFVReal
  []
  [permyy]
    type = MooseVariableFVReal
  []
  [permyz]
    type = MooseVariableFVReal
  []
  [permzx]
    type = MooseVariableFVReal
  []
  [permzy]
    type = MooseVariableFVReal
  []
  [permzz]
    type = MooseVariableFVReal
  []
  [poromat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permxxmat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permxymat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permxzmat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permyxmat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permyymat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permyzmat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permzxmat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permzymat]
    family = MONOMIAL
    order = CONSTANT
  []
  [permzzmat]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [poromat]
    type = ADPorousFlowPropertyAux
    property = porosity
    variable = poromat
  []
  [permxxmat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permxxmat
    column = 0
    row = 0
  []
  [permxymat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permxymat
    column = 1
    row = 0
  []
  [permxzmat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permxzmat
    column = 2
    row = 0
  []
  [permyxmat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permyxmat
    column = 0
    row = 1
  []
  [permyymat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permyymat
    column = 1
    row = 1
  []
  [permyzmat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permyzmat
    column = 2
    row = 1
  []
  [permzxmat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permzxmat
    column = 0
    row = 2
  []
  [permzymat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permzymat
    column = 1
    row = 2
  []
  [permzzmat]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = permzzmat
    column = 2
    row = 2
  []
[]

[ICs]
  [poro]
    type = RandomIC
    seed = 0
    variable = poro
    max = 0.5
    min = 0.1
  []
  [permx]
    type = FunctionIC
    function = permx
    variable = permxx
  []
  [permy]
    type = FunctionIC
    function = permy
    variable = permyy
  []
  [permz]
    type = FunctionIC
    function = permz
    variable = permzz
  []
[]

[Functions]
  [permx]
    type = ParsedFunction
    expression = '(1+x)*1e-11'
  []
  [permy]
    type = ParsedFunction
    expression = '(1+y)*1e-11'
  []
  [permz]
    type = ParsedFunction
    expression = '(1+z)*1e-11'
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    variable = ppwater
  []
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    variable = ppwater
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'ppwater'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 1000
    viscosity = 1e-3
    thermal_expansion = 0
    cv = 2
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
  []
  [ppss]
    type = ADPorousFlow1PhaseFullySaturated
    porepressure = ppwater
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = poro
  []
  [permeability]
    type = ADPorousFlowPermeabilityConstFromVar
    perm_xx = permxx
    perm_yy = permyy
    perm_zz = permzz
  []
  [relperm_water]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]

[Postprocessors]
  [mass_ph0]
    type = FVPorousFlowFluidMass
    fluid_component = 0
    execute_on = 'initial timestep_end'
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 100
  dt = 100
[]

[Outputs]
  execute_on = 'initial timestep_end'
  exodus = true
  perf_graph = true
[]

(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_2phasePS_fv.i)

# Pressure pulse in 1D with 2 phases, 2components - transient using FV

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 10
    xmin = 0
    xmax = 100
  []
[]

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

[Variables]
  [ppwater]
    type = MooseVariableFVReal
    initial_condition = 2e6
  []
  [sgas]
    type = MooseVariableFVReal
    initial_condition = 0.3
  []
[]

[AuxVariables]
  [massfrac_ph0_sp0]
    type = MooseVariableFVReal
    initial_condition = 1
  []
  [massfrac_ph1_sp0]
    type = MooseVariableFVReal
    initial_condition = 0
  []
  [ppgas]
    type = MooseVariableFVReal
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = ppwater
  []
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    variable = ppwater
    fluid_component = 0
  []
  [mass1]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 1
    variable = sgas
  []
  [flux1]
    type = FVPorousFlowAdvectiveFlux
    variable = sgas
    fluid_component = 1
  []
[]

[AuxKernels]
  [ppgas]
    type = ADPorousFlowPropertyAux
    property = pressure
    phase = 1
    variable = ppgas
  []
[]

[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
    thermal_expansion = 0
    viscosity = 1e-3
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 2e7
    density0 = 1
    thermal_expansion = 0
    viscosity = 1e-5
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
  []
  [ppss]
    type = ADPorousFlow2PhasePS
    phase0_porepressure = ppwater
    phase1_saturation = sgas
    capillary_pressure = pc
  []
  [massfrac]
    type = ADPorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [simple_fluid0]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
  []
  [relperm_water]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
  [relperm_gas]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 1
    phase = 1
  []
[]

[FVBCs]
  [leftwater]
    type = FVDirichletBC
    boundary = left
    value = 3e6
    variable = ppwater
  []
  [rightwater]
    type = FVDirichletBC
    boundary = right
    value = 2e6
    variable = ppwater
  []
  [sgas]
    type = FVDirichletBC
    boundary = 'left right'
    value = 0.3
    variable = sgas
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1e3
  end_time = 1e4
[]

[VectorPostprocessors]
  [pp]
    type = ElementValueSampler
    sort_by = x
    variable = 'ppwater ppgas'
  []
[]

[Outputs]
  file_base = pressure_pulse_1d_2phasePS_fv
  print_linear_residuals = false
  [csv]
    type = CSV
    execute_on = final
  []
  exodus = true
[]

(modules/porous_flow/test/tests/gravity/grav02e_fv.i)

# Checking that gravity head is established in the transient situation when 0<=saturation<=1 (note the less-than-or-equal-to).
# 2phase (PS), 2components, constant capillary pressure, constant fluid bulk-moduli for each phase, constant viscosity,
# constant permeability, Corey relative permeabilities with no residual saturation

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmax = 100
[]

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

[Variables]
  [ppwater]
    type = MooseVariableFVReal
    initial_condition = 1.5e6
  []
  [sgas]
    type = MooseVariableFVReal

    initial_condition = 0.3
  []
[]

[AuxVariables]
  [massfrac_ph0_sp0]
    type = MooseVariableFVReal
    initial_condition = 1
  []
  [massfrac_ph1_sp0]
    type = MooseVariableFVReal
    initial_condition = 0
  []
  [ppgas]
    type = MooseVariableFVReal
  []
  [swater]
    type = MooseVariableFVReal
  []
  [relpermwater]
    type = MooseVariableFVReal
  []
  [relpermgas]
    type = MooseVariableFVReal
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = ppwater
  []
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 0
    variable = ppwater
  []
  [mass1]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 1
    variable = sgas
  []
  [flux1]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 1
    variable = sgas
  []
[]

[AuxKernels]
  [ppgas]
    type = ADPorousFlowPropertyAux
    property = pressure
    phase = 1
    variable = ppgas
    execute_on = 'initial timestep_end'
  []
  [swater]
    type = ADPorousFlowPropertyAux
    property = saturation
    phase = 0
    variable = swater
    execute_on = 'initial timestep_end'
  []
  [relpermwater]
    type = ADPorousFlowPropertyAux
    property = relperm
    phase = 0
    variable = relpermwater
    execute_on = 'initial timestep_end'
  []
  [relpermgas]
    type = ADPorousFlowPropertyAux
    property = relperm
    phase = 1
    variable = relpermgas
    execute_on = 'initial timestep_end'
  []
[]

[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-5
    thermal_expansion = 0
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
  []
  [ppss]
    type = ADPorousFlow2PhasePS
    phase0_porepressure = ppwater
    phase1_saturation = sgas
    capillary_pressure = pc
  []
  [massfrac]
    type = ADPorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  []
  [simple_fluid0]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1e-11 0 0 0 1e-11 0  0 0 1e-11'
  []
  [relperm_water]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm_gas]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
[]

[VectorPostprocessors]
  [vars]
    type = ElementValueSampler
    variable = 'ppgas ppwater sgas swater'
    sort_by = x
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 5e3
  nl_abs_tol = 1e-12
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
  []
[]

[Outputs]
  execute_on = 'final'
  perf_graph = true
  csv = true
[]

(modules/porous_flow/test/tests/gravity/grav01a_fv.i)

# Checking that gravity head is established using FV
# 1phase, vanGenuchten, constant fluid-bulk, constant viscosity, constant permeability, Corey relative perm
# fully saturated
# For better agreement with the analytical solution (ana_pp), just increase nx

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmin = -1
  xmax = 0
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    type = MooseVariableFVReal
  []
[]

[ICs]
  [p]
    type = RandomIC
    variable = pp
    min = 0
    max = 1
  []
[]

[FVKernels]
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pp
    gravity = '-1 0 0'
  []
[]

[Functions]
  [ana_pp]
    type = ParsedFunction
    symbol_names = 'g B p0 rho0'
    symbol_values = '1 1.2 0 1'
    expression = '-B*log(exp(-p0/B)+g*rho0*x/B)' # expected pp at base
  []
[]

[FVBCs]
  [z]
    type = FVDirichletBC
    variable = pp
    boundary = right
    value = 0
  []
[]

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

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

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
  []
  [ppss]
    type = ADPorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1 0 0  0 2 0  0 0 3'
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
[]

[Postprocessors]
  [pp_base]
    type = PointValue
    variable = pp
    point = '-1 0 0'
  []
  [pp_analytical]
    type = FunctionValuePostprocessor
    function = ana_pp
    point = '-1 0 0'
  []
  [pp_00]
    type = PointValue
    variable = pp
    point = '0 0 0'
  []
  [pp_01]
    type = PointValue
    variable = pp
    point = '-0.1 0 0'
  []
  [pp_02]
    type = PointValue
    variable = pp
    point = '-0.2 0 0'
  []
  [pp_03]
    type = PointValue
    variable = pp
    point = '-0.3 0 0'
  []
  [pp_04]
    type = PointValue
    variable = pp
    point = '-0.4 0 0'
  []
  [pp_05]
    type = PointValue
    variable = pp
    point = '-0.5 0 0'
  []
  [pp_06]
    type = PointValue
    variable = pp
    point = '-0.6 0 0'
  []
  [pp_07]
    type = PointValue
    variable = pp
    point = '-0.7 0 0'
  []
  [pp_08]
    type = PointValue
    variable = pp
    point = '-0.8 0 0'
  []
  [pp_09]
    type = PointValue
    variable = pp
    point = '-0.9 0 0'
  []
  [pp_10]
    type = PointValue
    variable = pp
    point = '-1 0 0'
  []
[]

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

[Executioner]
  type = Steady
  solve_type = Newton
[]

[Outputs]
  [csv]
    type = CSV
  []
[]

(modules/porous_flow/test/tests/dispersion/disp01_fv.i)

# Test dispersive part of FVPorousFlowDispersiveFlux kernel by setting diffusion
# coefficients to zero. A pressure gradient is applied over the mesh to give a
# uniform velocity. Gravity is set to zero.
# Mass fraction is set to 1 on the left hand side and 0 on the right hand side.

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 20
  xmax = 10
  bias_x = 1.1
[]

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

[Variables]
  [pp]
    type = MooseVariableFVReal
  []
  [massfrac0]
    type = MooseVariableFVReal
  []
[]

[AuxVariables]
  [velocity]
    family = MONOMIAL
    order = FIRST
  []
[]

[AuxKernels]
  [velocity]
    type = ADPorousFlowDarcyVelocityComponent
    variable = velocity
    component = x
  []
[]

[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = pic
  []
  [massfrac0]
    type = ConstantIC
    variable = massfrac0
    value = 0
  []
[]

[Functions]
  [pic]
    type = ParsedFunction
    expression = '1.1e5-x*1e3'
  []
[]

[FVBCs]
  [xleft]
    type = FVDirichletBC
    value = 1
    variable = massfrac0
    boundary = left
  []
  [xright]
    type = FVDirichletBC
    value = 0
    variable = massfrac0
    boundary = right
  []
  [pright]
    type = FVDirichletBC
    variable = pp
    boundary = right
    value = 1e5
  []
  [pleft]
    type = FVDirichletBC
    variable = pp
    boundary = left
    value = 1.1e5
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [adv0]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pp
  []
  [diff0]
    type = FVPorousFlowDispersiveFlux
    variable = pp
    disp_trans = 0
    disp_long = 0.2
  []
  [mass1]
    type = FVPorousFlowMassTimeDerivative
    fluid_component = 1
    variable = massfrac0
  []
  [adv1]
    type = FVPorousFlowAdvectiveFlux
    fluid_component = 1
    variable = massfrac0
  []
  [diff1]
    type = FVPorousFlowDispersiveFlux
    fluid_component = 1
    variable = massfrac0
    disp_trans = 0
    disp_long = 0.2
  []
[]

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

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 1e9
    density0 = 1000
    viscosity = 0.001
    thermal_expansion = 0
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
  []
  [ppss]
    type = ADPorousFlow1PhaseFullySaturated
    porepressure = pp
  []
  [massfrac]
    type = ADPorousFlowMassFraction
    mass_fraction_vars = massfrac0
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [poro]
    type = ADPorousFlowPorosityConst
    porosity = 0.3
  []
  [diff]
    type = ADPorousFlowDiffusivityConst
    diffusion_coeff = '0 0'
    tortuosity = 0.1
  []
  [relp]
    type = ADPorousFlowRelativePermeabilityConst
    phase = 0
  []
  [permeability]
    type = ADPorousFlowPermeabilityConst
    permeability = '1e-9 0 0 0 1e-9 0 0 0 1e-9'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
    petsc_options_value = 'gmres      asm      lu           NONZERO'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 3e2
  dtmax = 100
  nl_abs_tol = 1e-12
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 2
    cutback_factor = 0.5
    dt = 10
  []
[]

[VectorPostprocessors]
  [xmass]
    type = ElementValueSampler
    sort_by = id
    variable = 'massfrac0 velocity'
  []
[]

[Outputs]
  [out]
    type = CSV
    execute_on = final
  []
[]

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

# Testing permeability from porosity
# Trivial test, checking calculated permeability is correct
# k = k_anisotropic * B * exp(A * phi)

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 3
    xmin = 0
    xmax = 3
  []
[]

[GlobalParams]
  block = 0
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    type = MooseVariableFVReal
    [FVInitialCondition]
      type = FVConstantIC
      value = 0
    []
  []
[]

[FVKernels]
  [flux]
    type = FVPorousFlowAdvectiveFlux
    gravity = '0 0 0'
    variable = pp
  []
[]

[FVBCs]
  [ptop]
    type = FVDirichletBC
    variable = pp
    boundary = right
    value = 0
  []
  [pbase]
    type = FVDirichletBC
    variable = pp
    boundary = left
    value = 1
  []
[]

[AuxVariables]
  [poro]
    type = MooseVariableFVReal
  []
  [perm_x]
    type = MooseVariableFVReal
  []
  [perm_y]
    type = MooseVariableFVReal
  []
  [perm_z]
    type = MooseVariableFVReal
  []
[]

[AuxKernels]
  [poro]
    type = ADPorousFlowPropertyAux
    property = porosity
    variable = poro
  []
  [perm_x]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_x
    row = 0
    column = 0
  []
  [perm_y]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_y
    row = 1
    column = 1
  []
  [perm_z]
    type = ADPorousFlowPropertyAux
    property = permeability
    variable = perm_z
    row = 2
    column = 2
  []
[]

[Postprocessors]
  [perm_x_bottom]
    type = PointValue
    variable = perm_x
    point = '0 0 0'
  []
  [perm_y_bottom]
    type = PointValue
    variable = perm_y
    point = '0 0 0'
  []
  [perm_z_bottom]
    type = PointValue
    variable = perm_z
    point = '0 0 0'
  []
  [perm_x_top]
    type = PointValue
    variable = perm_x
    point = '3 0 0'
  []
  [perm_y_top]
    type = PointValue
    variable = perm_y
    point = '3 0 0'
  []
  [perm_z_top]
    type = PointValue
    variable = perm_z
    point = '3 0 0'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    # unimportant in this fully-saturated test
    m = 0.8
    alpha = 1e-4
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2.2e9
    viscosity = 1e-3
    density0 = 1000
    thermal_expansion = 0
  []
[]

[Materials]
  [permeability]
    type = ADPorousFlowPermeabilityExponential
    k_anisotropy = '1 0 0  0 2 0  0 0 0.1'
    poroperm_function = exp_k
    A = 10
    B = 1e-8
  []
  [temperature]
    type = ADPorousFlowTemperature
  []
  [massfrac]
    type = ADPorousFlowMassFraction
  []
  [eff_fluid_pressure]
    type = ADPorousFlowEffectiveFluidPressure
  []
  [ppss]
    type = ADPorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [simple_fluid]
    type = ADPorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = 0.1
  []
  [relperm]
    type = ADPorousFlowRelativePermeabilityCorey
    n = 0 # unimportant in this fully-saturated situation
    phase = 0
  []
[]

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

[Executioner]
  type = Steady
  solve_type = Newton
  l_tol = 1E-5
  nl_abs_tol = 1E-3
  nl_rel_tol = 1E-8
  l_max_its = 200
  nl_max_its = 400
[]

[Outputs]
  file_base = 'PermFromPoro03_out'
  csv = true
  execute_on = 'timestep_end'
[]

(modules/porous_flow/examples/fluidflower/fluidflower.i)

# FluidFlower International Benchmark study model
# CSIRO 2023
#
# This example can be used to reproduce the results presented by the
# CSIRO team as part of this benchmark study. See
# Green, C., Jackson, S.J., Gunning, J., Wilkins, A. and Ennis-King, J.,
# 2023. Modelling the FluidFlower: Insights from Characterisation and
# Numerical Predictions. Transport in Porous Media.
#
# This example takes a long time to run! The large density contrast
# between the gas phase CO2 and the water makes convergence very hard,
# so small timesteps must be taken during injection.
#
# This example uses a simplified mesh in order to be run during the
# automated testing. To reproduce the results of the benchmark study,
# replace the simple layered input mesh with the one located in the
# large_media submodule.
#
# The mesh file contains:
# - porosity as given by FluidFlower description
# - permeability as given by FluidFlower description
# - subdomain ids for each sand type
#
# The nominal thickness of the FluidFlower tank is 19mm. To keep masses consistent
# with the experiment, porosity and permeability are multiplied by the thickness
thickness = 0.019
#
# Properties associated with each sand type associated with mesh block ids
#
# block 0 - ESF (very fine sand)
sandESF = '0 10 20'
sandESF_pe = 1471.5
sandESF_krg = 0.09
sandESF_swi = 0.32
sandESF_krw = 0.71
sandESF_sgi = 0.14

# block 1 - C - Coarse lower
sandC = '1 21'
sandC_pe = 294.3
sandC_krg = 0.05
sandC_swi = 0.14
sandC_krw = 0.93
sandC_sgi = 0.1

# block 2 - D - Coarse upper
sandD = '2 22'
sandD_pe = 98.1
sandD_krg = 0.02
sandD_swi = 0.12
sandD_krw = 0.95
sandD_sgi = 0.08

# block 3 - E - Very Coarse lower
sandE = '3 13 23'
sandE_pe = 10
sandE_krg = 0.1
sandE_swi = 0.12
sandE_krw = 0.93
sandE_sgi = 0.06

# block 4 - F - Very Coarse upper
sandF = '4 14 24 34'
sandF_pe = 10
sandF_krg = 0.11
sandF_swi = 0.12
sandF_krw = 0.72
sandF_sgi = 0.13

# block 5 - G - Flush Zone
sandG = '5 15 35'
sandG_pe = 10
sandG_krg = 0.16
sandG_swi = 0.1
sandG_krw = 0.75
sandG_sgi = 0.06

# block 6 - Fault 1 - Heterogeneous
fault1 = '6 26'
fault1_pe = 10
fault1_krg = 0.16
fault1_swi = 0.1
fault1_krw = 0.75
fault1_sgi = 0.06

# block 7 - Fault 2 - Impermeable
# Note: this fault has been removed from the mesh (no elements in this region)

# block 8 - Fault 3 - Homogeneous
fault3 = '8'
fault3_pe = 10
fault3_krg = 0.16
fault3_swi = 0.1
fault3_krw = 0.75
fault3_sgi = 0.06

# Top layer
top_layer = '9'

# Boxes A, B an C used to report values (sg, sgr, xco2, etc)
boxA = '10 13 14 15 34 35'
boxB = '20 21 22 23 24 26'
boxC = '34 35'

# Furthermore, the seal sand unit in boxes A and B
seal_boxA = '10'
seal_boxB = '20'

# CO2 injection details:
# CO2 density ~1.8389 kg/m3 at 293.15 K, 1.01325e5 Pa
# Injection in Port (9, 3) for 5 hours.
# Injection in Port (17, 7) for 2:45 hours.
# Injection of 10 ml/min = 0.1666 ml/s = 1.666e-7 m3/s = ~3.06e-7 kg/s.
# Total mass of CO2 injected ~ 8.5g.
inj_rate = 3.06e-7

[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = 'fluidflower_test.e'
    # file = '../../../../large_media/porous_flow/examples/fluidflower/fluidflower.e'
    use_for_exodus_restart = true
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 -9.81 0'
  temperature = temperature
  log_extension = false
[]

[Variables]
  [pgas]
    family = MONOMIAL
    order = CONSTANT
    fv = true
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    scaling = 1e4
  []
[]

[AuxVariables]
  [xnacl]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 0.0055
  []
  [temperature]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 20
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_from_file_var = porosity
  []
  [porosity_times_thickness]
    family = MONOMIAL
    order = CONSTANT
    fv = true
  []
  [permeability]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_from_file_var = permeability
  []
  [permeability_times_thickness]
    family = MONOMIAL
    order = CONSTANT
    fv = true
  []
  [saturation_water]
    family = MONOMIAL
    order = CONSTANT
  []
  [saturation_gas]
    family = MONOMIAL
    order = CONSTANT
  []
  [pressure_water]
    family = MONOMIAL
    order = CONSTANT
  []
  [pc]
    family = MONOMIAL
    order = CONSTANT
  []
  [x0_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x0_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [porosity_times_thickness]
    type = ParsedAux
    variable = porosity_times_thickness
    coupled_variables = porosity
    expression = 'porosity * ${thickness}'
    execute_on = 'initial'
  []
  [permeability_times_thickness]
    type = ParsedAux
    variable = permeability_times_thickness
    coupled_variables = permeability
    expression = 'permeability * ${thickness}'
    execute_on = 'initial'
  []
  [pressure_water]
    type = ADPorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_water]
    type = ADPorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_gas]
    type = ADPorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [density_water]
    type = ADPorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = ADPorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [x1_water]
    type = ADPorousFlowPropertyAux
    variable = x1_water
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = 'initial timestep_end'
  []
  [x1_gas]
    type = ADPorousFlowPropertyAux
    variable = x1_gas
    property = mass_fraction
    phase = 1
    fluid_component = 1
    execute_on = 'initial timestep_end'
  []
  [x0_water]
    type = ADPorousFlowPropertyAux
    variable = x0_water
    property = mass_fraction
    phase = 0
    fluid_component = 0
    execute_on = 'initial timestep_end'
  []
  [x0_gas]
    type = ADPorousFlowPropertyAux
    variable = x0_gas
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = 'initial timestep_end'
  []
  [pc]
    type = ADPorousFlowPropertyAux
    variable = pc
    property = capillary_pressure
    execute_on = 'initial timestep_end'
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [diff0]
    type = FVPorousFlowDispersiveFlux
    variable = pgas
    fluid_component = 0
    disp_long = '0 0'
    disp_trans = '0 0'
  []
  [mass1]
    type = FVPorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [flux1]
    type = FVPorousFlowAdvectiveFlux
    variable = z
    fluid_component = 1
  []
  [diff1]
    type = FVPorousFlowDispersiveFlux
    variable = z
    fluid_component = 1
    disp_long = '0 0'
    disp_trans = '0 0'
  []
[]

[DiracKernels]
  [injector1]
    type = ConstantPointSource
    point = '0.9 0.3 0'
    value = ${inj_rate}
    variable = z
  []
  [injector2]
    type = ConstantPointSource
    point = '1.7 0.7 0'
    value = ${inj_rate}
    variable = z
  []
[]

[Controls]
  [injection1]
    type = ConditionalFunctionEnableControl
    enable_objects = 'DiracKernels::injector1'
    conditional_function = injection_schedule1
  []
  [injection2]
    type = ConditionalFunctionEnableControl
    enable_objects = 'DiracKernels::injector2'
    conditional_function = injection_schedule2
  []
[]

[Functions]
  [initial_p]
    type = ParsedFunction
    symbol_names = 'p0 g H rho0'
    symbol_values = '101.325e3 9.81 1.5 1002'
    expression = 'p0 + rho0 * g * (H - y)'
  []
  [injection_schedule1]
    type = ParsedFunction
    expression = 'if(t >= 0 & t <= 1.8e4, 1, 0)'
  []
  [injection_schedule2]
    type = ParsedFunction
    expression = 'if(t >= 8.1e3 & t <= 1.8e4, 1, 0)'
  []
[]

[ICs]
  [p]
    type = FunctionIC
    variable = pgas
    function = initial_p
  []
[]

[FVBCs]
  [pressure_top]
    type = FVPorousFlowAdvectiveFluxBC
    boundary = top
    porepressure_value = 1.01325e5
    variable = pgas
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
  [watertab]
    type = TabulatedBicubicFluidProperties
    fp = water
    save_file = false
    pressure_min = 1e5
    pressure_max = 1e6
    temperature_min = 290
    temperature_max = 300
    num_p = 20
    num_T = 10
  []
  [co2]
    type = CO2FluidProperties
  []
  [co2tab]
    type = TabulatedBicubicFluidProperties
    fp = co2
    save_file = false
    pressure_min = 1e5
    pressure_max = 1e6
    temperature_min = 290
    temperature_max = 300
    num_p = 20
    num_T = 10
  []
  [brine]
    type = BrineFluidProperties
    water_fp = watertab
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [sandESF_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandESF_pe}
    lambda = 2
    block = ${sandESF}
    pc_max = 1e4
    sat_lr = ${sandESF_swi}
  []
  [sandC_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandC_pe}
    lambda = 2
    block = ${sandC}
    pc_max = 1e4
    sat_lr = ${sandC_swi}
  []
  [sandD_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandD_pe}
    lambda = 2
    block = ${sandD}
    pc_max = 1e4
    sat_lr = ${sandD_swi}
  []
  [sandE_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandE_pe}
    lambda = 2
    block = ${sandE}
    pc_max = 1e4
    sat_lr = ${sandE_swi}
  []
  [sandF_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandF_pe}
    lambda = 2
    block = ${sandF}
    pc_max = 1e4
    sat_lr = ${sandF_swi}
  []
  [sandG_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandG_pe}
    lambda = 2
    block = ${sandG}
    pc_max = 1e4
    sat_lr = ${sandG_swi}
  []
  [fault1_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${fault1_pe}
    lambda = 2
    block = ${fault1}
    pc_max = 1e4
    sat_lr = ${fault1_swi}
  []
  [fault3_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${fault3_pe}
    lambda = 2
    block = ${fault3}
    pc_max = 1e4
    sat_lr = ${fault3_swi}
  []
  [top_layer_pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
    block =  ${top_layer}
  []
  [sandESF_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandESF_pc
  []
  [sandC_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandC_pc
  []
  [sandD_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandD_pc
  []
  [sandE_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandE_pc
  []
  [sandF_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandF_pc
  []
  [sandG_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandG_pc
  []
  [fault1_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = fault1_pc
  []
  [fault3_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = fault3_pc
  []
  [top_layer_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = top_layer_pc
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
    temperature = temperature
  []
  [sandESF_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandESF_fs
    capillary_pressure = sandESF_pc
    block = ${sandESF}
  []
  [sandC_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandC_fs
    capillary_pressure = sandC_pc
    block = ${sandC}
  []
  [sandD_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandD_fs
    capillary_pressure = sandD_pc
    block = ${sandD}
  []
  [sandE_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandE_fs
    capillary_pressure = sandE_pc
    block = ${sandE}
  []
  [sandF_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandF_fs
    capillary_pressure = sandF_pc
    block = ${sandF}
  []
  [sandG_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandG_fs
    capillary_pressure = sandG_pc
    block = ${sandG}
  []
  [fault1_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = fault1_fs
    capillary_pressure = fault1_pc
    block = ${fault1}
  []
  [fault3_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = fault3_fs
    capillary_pressure = fault3_pc
    block = ${fault3}
  []
  [top_layer_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = top_layer_fs
    capillary_pressure = top_layer_pc
    block = ${top_layer}
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = porosity_times_thickness
  []
  [permeability]
    type = ADPorousFlowPermeabilityConstFromVar
    perm_xx = permeability_times_thickness
    perm_yy = permeability_times_thickness
    perm_zz = permeability_times_thickness
  []
  [diffcoeff]
    type = ADPorousFlowDiffusivityConst
    tortuosity = '1 1'
    diffusion_coeff = '2e-9 2e-9 0 0'
  []
  [sandESF_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandESF_swi}
    sum_s_res = ${fparse sandESF_sgi + sandESF_swi}
    scaling = ${sandESF_krw}
    block = ${sandESF}
  []
  [sandESF_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandESF_sgi}
    sum_s_res = ${fparse sandESF_sgi + sandESF_swi}
    scaling = ${sandESF_krg}
    block = ${sandESF}
  []
  [sandC_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandC_swi}
    sum_s_res = ${fparse sandC_sgi + sandC_swi}
    scaling = ${sandC_krw}
    block = ${sandC}
  []
  [sandC_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandC_sgi}
    sum_s_res = ${fparse sandC_sgi + sandC_swi}
    scaling = ${sandC_krg}
    block = ${sandC}
  []
  [sandD_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandD_swi}
    sum_s_res = ${fparse sandD_sgi + sandD_swi}
    scaling = ${sandD_krw}
    block = ${sandD}
  []
  [sandD_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandD_sgi}
    sum_s_res = ${fparse sandD_sgi + sandD_swi}
    scaling = ${sandD_krg}
    block = ${sandD}
  []
  [sandE_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandE_swi}
    sum_s_res = ${fparse sandE_sgi + sandE_swi}
    scaling = ${sandE_krw}
    block = ${sandE}
  []
  [sandE_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandE_sgi}
    sum_s_res = ${fparse sandE_sgi + sandE_swi}
    scaling = ${sandE_krg}
    block = ${sandE}
  []
  [sandF_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandF_swi}
    sum_s_res = ${fparse sandF_sgi + sandF_swi}
    scaling = ${sandF_krw}
    block = ${sandF}
  []
  [sandF_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandF_sgi}
    sum_s_res = ${fparse sandF_sgi + sandF_swi}
    scaling = ${sandF_krg}
    block = ${sandF}
  []
  [sandG_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandG_swi}
    sum_s_res = ${fparse sandG_sgi + sandG_swi}
    scaling = ${sandG_krw}
    block = ${sandG}
  []
  [sandG_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandG_sgi}
    sum_s_res = ${fparse sandG_sgi + sandG_swi}
    scaling = ${sandG_krg}
    block = ${sandG}
  []
  [fault1_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${fault1_swi}
    sum_s_res = ${fparse fault1_sgi + fault1_swi}
    scaling = ${fault1_krw}
    block = ${fault1}
  []
  [fault1_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${fault1_sgi}
    sum_s_res = ${fparse fault1_sgi + fault1_swi}
    scaling = ${fault1_krg}
    block = ${fault1}
  []
  [fault3_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${fault3_swi}
    sum_s_res = ${fparse fault3_sgi + fault3_swi}
    scaling = ${fault3_krw}
    block = ${fault3}
  []
  [fault3_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${fault3_sgi}
    sum_s_res = ${fparse fault3_sgi + fault3_swi}
    scaling = ${fault3_krg}
    block = ${fault3}
  []
  [top_layer_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    block = ${top_layer}
  []
  [top_layer_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    block = ${top_layer}
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options = '-ksp_snes_ew'
    petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package -sub_pc_factor_shift_type'
    petsc_options_value = 'gmres lu mumps NONZERO'
    # petsc_options_iname = '-ksp_type -pc_type -pc_hypre_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels -ksp_gmres_restart'
    # petsc_options_value = 'gmres hypre boomeramg lu NONZERO 4 301'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dtmax = 60
  start_time = 0
  end_time = 4.32e5
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8
  nl_max_its = 15
  l_tol = 1e-5
  l_abs_tol = 1e-8
  # line_search = none # Can be a useful option for this problem
  [TimeSteppers]
    [time]
      type = FunctionDT
      growth_factor = 2
      cutback_factor_at_failure = 0.5
      function = 'if(t<1.8e4, 2, if(t<3.6e4, 20, 60))'
    []
  []
[]

[Postprocessors]
  [p_5_3]
    type = PointValue
    variable = pgas
    point = '0.5 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_5_3_w]
    type = PointValue
    variable = pressure_water
    point = '0.5 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_5_7]
    type = PointValue
    variable = pgas
    point = '0.5 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_5_7_w]
    type = PointValue
    variable = pressure_water
    point = '0.5 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_9_3]
    type = PointValue
    variable = pgas
    point = '0.9 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_9_3_w]
    type = PointValue
    variable = pressure_water
    point = '0.9 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_15_5]
    type = PointValue
    variable = pgas
    point = '1.5 0.5 0'
    execute_on = 'initial timestep_end'
  []
  [p_15_5_w]
    type = PointValue
    variable = pressure_water
    point = '1.5 0.5 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_7]
    type = PointValue
    variable = pgas
    point = '1.7 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_7_w]
    type = PointValue
    variable = pressure_water
    point = '1.7 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_11]
    type = PointValue
    variable = pgas
    point = '1.7 1.1 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_11_w]
    type = PointValue
    variable = pressure_water
    point = '1.7 1.1 0'
    execute_on = 'initial timestep_end'
  []
  [x0mass]
    type = FVPorousFlowFluidMass
    fluid_component = 0
    phase = '0 1'
    execute_on = 'initial timestep_end'
  []
  [x1mass]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    execute_on = 'initial timestep_end'
  []
  [x1gas]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '1'
    execute_on = 'initial timestep_end'
  []
  [boxA]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${boxA}
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandESF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandESF_sgi}
    block = 10
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandE]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandE_sgi}
    block = 13
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandF_sgi}
    block = '14 34'
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandG]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandG_sgi}
    block = '15 35'
    execute_on = 'initial timestep_end'
  []
  [imm_A]
    type = LinearCombinationPostprocessor
    pp_names = 'imm_A_sandESF imm_A_sandE imm_A_sandF imm_A_sandG'
    pp_coefs = '1 1 1 1'
    execute_on = 'initial timestep_end'
  []
  [diss_A]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 0
    block = ${boxA}
    execute_on = 'initial timestep_end'
  []
  [seal_A]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${seal_boxA}
    execute_on = 'initial timestep_end'
  []
  [boxB]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${boxB}
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandESF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandESF_sgi}
    block = 20
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandC]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandC_sgi}
    block = 21
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandD]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandD_sgi}
    block = 22
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandE]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandE_sgi}
    block = 23
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandF_sgi}
    block = 24
    execute_on = 'initial timestep_end'
  []
  [imm_B_fault1]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${fault1_sgi}
    block = 26
    execute_on = 'initial timestep_end'
  []
  [imm_B]
    type = LinearCombinationPostprocessor
    pp_names = 'imm_B_sandESF imm_B_sandC imm_B_sandD imm_B_sandE imm_B_sandF imm_B_fault1'
    pp_coefs = '1 1 1 1 1 1'
    execute_on = 'initial timestep_end'
  []
  [diss_B]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 0
    block = ${boxB}
    execute_on = 'initial timestep_end'
  []
  [seal_B]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${seal_boxB}
    execute_on = 'initial timestep_end'
  []
  [boxC]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0'
    block = ${boxC}
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  # exodus = true
  [csv]
     type = CSV
  []
[]

Input Parameters
Input Files