How to use Kuzmin-Turek stabilization in PorousFlow simulations

This page is part of a set of pages devoted to discussions of numerical stabilization in PorousFlow. See:

Kuzmin and Turek (Kuzmin and Turek, 2004) describe a method of stabilising advection while minimising artificial numerical diffusion. In this page "Kuzmin and Turek" is abbreviated to "KT". Basic features of KT's approach are detailed in a worked example and the results are compared with full upwinding, RDG and no stabilization in the numerical diffusion page.

This page is for users who want to use KT stabilization in PorousFlow simulations.

note

The Kuzmin-Turek stabilization is new and users are strongly encouraged to experiment with it and report their findings to the moose-users gmail group to iron out any problems and so we can collectively gain experience.

Using the Action system

Both the PorousFlowFullySaturated and PorousFlowUnsaturated Actions support KT stabilization (the PorousFlowBasicTHM Action does not support any numerical stabilization as it is typically unneeded). KT stabilization may be included simply by specifying the stabilization and flux_limiter_type parameters. For instance:

[PorousFlowFullySaturated<<<{"href": "../../syntax/PorousFlowFullySaturated/index.html"}>>>]
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  coupling_type<<<{"description": "The type of simulation.  For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient.  For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = Hydro
  gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = the_simple_fluid
  mass_fraction_vars<<<{"description": "List of variables that represent the mass fractions.  With only one fluid component, this may be left empty.  With N fluid components, the format is 'f_0 f_1 f_2 ... f_(N-1)'.  That is, the N^th component need not be specified because f_N = 1 - (f_0 + f_1 + ... + f_(N-1)).  It is best numerically to choose the N-1 mass fraction variables so that they represent the fluid components with small concentrations.  This Action will associated the i^th mass fraction variable to the equation for the i^th fluid component, and the pressure variable to the N^th fluid component."}>>> = tracer_concentration
  stabilization<<<{"description": "Numerical stabilization used.  'Full' means full upwinding.  'KT' means FEM-TVD stabilization of Kuzmin-Turek"}>>> = KT
  flux_limiter_type<<<{"description": "Type of flux limiter to use if stabilization=KT.  'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = superbee
[]

and

[PorousFlowUnsaturated<<<{"href": "../../syntax/PorousFlowUnsaturated/index.html"}>>>]
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  coupling_type<<<{"description": "The type of simulation.  For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient.  For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = Hydro
  gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = the_simple_fluid
  relative_permeability_exponent<<<{"description": "Relative permeability exponent"}>>> = 3
  relative_permeability_type<<<{"description": "Type of relative-permeability function.  FLAC relperm = (1+m)S^m - mS^(1+m).  Corey relperm = S^m.  m is the exponent.  Here S = (saturation - residual)/(1 - residual)"}>>> = Corey
  residual_saturation<<<{"description": "Residual saturation to use in the relative permeability expression"}>>> = 0.1
  van_genuchten_alpha<<<{"description": "Van Genuchten alpha parameter used to determine saturation from porepressure"}>>> = 1E-6
  van_genuchten_m<<<{"description": "Van Genuchten m parameter used to determine saturation from porepressure"}>>> = 0.6
  stabilization<<<{"description": "Numerical stabilization used.  'Full' means full upwinding.  'KT' means FEM-TVD stabilization of Kuzmin-Turek"}>>> = KT
  flux_limiter_type<<<{"description": "Type of flux limiter to use if stabilization=KT.  'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = None
[]

Using Kernels and UserObjects

It is also fairly straightforward to swap from full-upwinding (or no stabilization) to KT stabilization. Only two changes are needed in the input file.

Kernels

Instead of PorousFlowAdvectiveFlux or PorousFlowHeatAdvection or another type of PorousFlow Kernel, the PorousFlowFluxLimitedTVDAdvection Kernel must be used. For instance, a full-upwind simulation's Kernels are:

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [mass_dot]
    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
  []
  [advection]
    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
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
  []
  [heat_advection]
    type = PorousFlowHeatAdvection<<<{"description": "Fully-upwinded heat flux, advected by the fluid", "href": "../../source/kernels/PorousFlowHeatAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
  []
[]

while an identical simulation using KT stabilization has Kernels:

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [mass_dot]
    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
  []
  [fluid_advection]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = fluid_advective_flux
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
  []
  [heat_advection]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = heat_advective_flux
  []
[]

Notice that the PorousFlowFluxLimitedTVDAdvection Kernels require an advective_flux_calculator input. This is a UserObject, which is the subject of the next section.

UserObjects

The set of PorousFlow advective flux calculator UserObjects compute the Residual and Jacobian entries that are used by the PorousFlowFluxLimitedTVDAdvection Kernels. The reason that these entries are computed by UserObjects instead of Kernels is that the KT scheme requires information about porepressures, etc, that are 2 nodes away from the current node. This information gets assembled and processed to determine upwinding directions and suitable stabilization approaches before the residual and Jacobian are computed. The process is explained in the worked example of KT stabilization page. The upshot is that a UserObject must be employed instead of the usual MOOSE Kernel-only approach.

Table 1: Available types of PorousFlow advective flux UserObjects

Equation	UserObject	Use case
$var element = document.getElementById("moose-equation-b2122058-c605-4100-bf2e-d72553a7125d");katex.render("-\\frac{\\rho}{\\mu}k_{ij}(\\nabla_{j}P - \\rho g_{j})", element, {displayMode:false,throwOnError:false});$	PorousFlowAdvectiveFluxCalculatorSaturated	1-phase, 1-component, fully-saturated
$var element = document.getElementById("moose-equation-b02888c1-cb35-418f-b5d1-688c403b101f");katex.render("-\\chi^{\\kappa}\\frac{\\rho}{\\mu}k_{ij}(\\nabla_{j}P - \\rho g_{j})", element, {displayMode:false,throwOnError:false});$	PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent	1-phase, multi-component, fully-saturated
$var element = document.getElementById("moose-equation-b2fc8ada-6236-4c2c-9fcf-278dba497da4");katex.render("-k_{\\mathrm{r}}\\frac{\\rho}{\\mu}k_{ij}(\\nabla_{j}P - \\rho g_{j})", element, {displayMode:false,throwOnError:false});$	PorousFlowAdvectiveFluxCalculatorUnsaturated	1-phase, 1-component, unsaturated
$var element = document.getElementById("moose-equation-aacf0b92-1b24-45e7-930a-c493f31ca85f");katex.render("-k_{\\mathrm{r,}\\beta}\\frac{\\rho_{\\beta}}{\\mu_{\\beta}}k_{ij}(\\nabla_{j}P_{\\beta} - \\rho_{\\beta} g_{j})", element, {displayMode:false,throwOnError:false});$	PorousFlowAdvectiveFluxCalculatorUnsaturated	multi-phase, multi-component, immiscible
$var element = document.getElementById("moose-equation-7c997317-f8ce-448e-a8b2-6cdf781700ab");katex.render("-k_{\\mathrm{r,}\\beta}\\chi_{\\beta}^{\\kappa}\\frac{\\rho_{\\beta}}{\\mu_{\\beta}}k_{ij}(\\nabla_{j}P_{\\beta} - \\rho_{\\beta} g_{j})", element, {displayMode:false,throwOnError:false});$	PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent	multi-phase, multi-component
$var element = document.getElementById("moose-equation-42a4cc29-4094-4905-8ef5-72b1dcccd11a");katex.render("-h\\frac{\\rho}{\\mu}k_{ij}(\\nabla_{j}P - \\rho g_{j})", element, {displayMode:false,throwOnError:false});$	PorousFlowAdvectiveFluxCalculatorSaturatedHeat	heat, 1-phase
$var element = document.getElementById("moose-equation-2e0686f4-93a3-4284-9f23-fcec86b379d5");katex.render("-h_{\\beta}k_{\\mathrm{r,}\\beta}\\frac{\\rho_{\\beta}}{\\mu_{\\beta}}k_{ij}(\\nabla_{j}P_{\\beta} - \\rho_{\\beta} g_{j})", element, {displayMode:false,throwOnError:false});$	PorousFlowAdvectiveFluxCalculatorUnsaturatedHeat	heat, multi-phase

The notation is detailed in the nomenclature: briefly, a $var element = document.getElementById("moose-equation-9442da44-b52e-453d-96e9-e5572cc908c7");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ subscript denotes a phase, a $var element = document.getElementById("moose-equation-7bd6d7fa-cce5-4a84-8131-71267fd96a61");katex.render("\\kappa", element, {displayMode:false,throwOnError:false});$ superscript denotes a fluid component, $var element = document.getElementById("moose-equation-b889123f-4124-4e47-bda0-073dec0abdc3");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is a fluid density, $var element = document.getElementById("moose-equation-6c8dd278-7294-42d0-95db-d8a18d3ed9d2");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ a fluid viscosity, $var element = document.getElementById("moose-equation-41c1c341-3eac-4e34-a706-6cb1004a7ea3");katex.render("k_{ij}", element, {displayMode:false,throwOnError:false});$ the permeability tensor, $var element = document.getElementById("moose-equation-7b09d40f-01fb-4ac2-b023-8ffb79ab11be");katex.render("P", element, {displayMode:false,throwOnError:false});$ a fluid porepressure, $var element = document.getElementById("moose-equation-46a51bee-0942-4e7c-8874-e442b1cd1dc7");katex.render("g_{j}", element, {displayMode:false,throwOnError:false});$ the acceleration of gravity, $var element = document.getElementById("moose-equation-e5f6137c-7111-42e0-8a9b-d5ab872d1ce4");katex.render("\\chi", element, {displayMode:false,throwOnError:false});$ a mass fraction, $var element = document.getElementById("moose-equation-5cc25cb9-ae93-484f-88f6-34e9851dbcf5");katex.render("k_{\\mathrm{r}}", element, {displayMode:false,throwOnError:false});$ a relative permeability, and $var element = document.getElementById("moose-equation-faf0704d-b0e2-42de-a6f1-dcf889908ac8");katex.render("h", element, {displayMode:false,throwOnError:false});$ an enthalpy.

Which UserObjects you need depends on your simulation. For example, in the case of a single-phase, single-component, fully-saturated anisothermal simulation we need two UserObjects: one for the advection of the fluid, and one for the advection of the heat.

  [fluid_advective_flux]
    type = PorousFlowAdvectiveFluxCalculatorSaturated
    flux_limiter_type = superbee
  []
  [heat_advective_flux]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedHeat
    flux_limiter_type = superbee
  []
[]

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

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1.0
    density = 125
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 2 0 0 0 3'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [PS]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres bjacobi 1E-15 1E-10 10000'
  []
[]

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

[VectorPostprocessors]
  [T]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '1 0 0'
    num_points = 51
    sort_by = x
    variable = temp
  []
[]

[Outputs]
  file_base = heat_advection_1d_KT
  [csv]
    type = CSV
    sync_times = '0.1 0.6'
    sync_only = true
  []
[]

Clicking on the links in the above table will provide you with more examples. For instance the following set of Kernels and UserObjects illustrates the case of a single tracer in a fully-saturated single-phase fluid:

[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"}>>> = tracer
  []
  [flux0]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tracer
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = advective_flux_calculator_0
  []
  [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"}>>> = porepressure
  []
  [flux1]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = porepressure
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = advective_flux_calculator_1
  []
[]

  [advective_flux_calculator_0]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
    flux_limiter_type = superbee
    fluid_component = 0
  []
  [advective_flux_calculator_1]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
    flux_limiter_type = superbee
    fluid_component = 1
  []
[]

note

For multi-phase situations you will need an advective flux calculator for each phase and each component, unless the phases are immiscible (each component exists in one phase only).

For example, in the case of 2 phases with 2 components, each potentially existing in both phases, there are 2 PorousFlow governing equations: one for each component. The equation for fluid component zero contains contributions from both phase 0 and phase 1. The equation for fluid component one contains contributions from both phase 0 and phase 1. So the Kernels will look like

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [mass_component0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
  []
  [flux_component0_phase0]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = afc_component0_phase0
  []
  [flux_component0_phase1]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = afc_component0_phase1
  []
  [mass_component1]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = sgas
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 1
  []
  [flux_component1_phase0]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = sgas
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = afc_component1_phase0
  []
  [flux_component1_phase1]
    type = PorousFlowFluxLimitedTVDAdvection<<<{"description": "Advective flux of fluid species or heat using the Flux Limited TVD scheme invented by Kuzmin and Turek", "href": "../../source/kernels/PorousFlowFluxLimitedTVDAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = sgas
    advective_flux_calculator<<<{"description": "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection describes a movement of a fluid component in a fluid phase, or movement of heat energy in a fluid phase"}>>> = afc_component1_phase1
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [ppgas]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = pressure
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 1
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = ppgas
  []
[]

[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
  []
  [afc_component0_phase0]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent<<<{"description": "Computes the advective flux of fluid of given phase and component.  Hence this UserObject is relevant to multi-phase, multi-component situations.  Explicitly, the UserObject computes (mass_fraction * density * relative_permeability / viscosity) * (- permeability * (grad(P) - density * gravity)), using the Kuzmin-Turek FEM-TVD multidimensional stabilization scheme", "href": "../../source/userobjects/PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this UserObject"}>>> = 0
    phase<<<{"description": "The index corresponding to the phase for this UserObject"}>>> = 0
    flux_limiter_type<<<{"description": "Type of flux limiter to use.  'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = superbee
  []
  [afc_component0_phase1]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent<<<{"description": "Computes the advective flux of fluid of given phase and component.  Hence this UserObject is relevant to multi-phase, multi-component situations.  Explicitly, the UserObject computes (mass_fraction * density * relative_permeability / viscosity) * (- permeability * (grad(P) - density * gravity)), using the Kuzmin-Turek FEM-TVD multidimensional stabilization scheme", "href": "../../source/userobjects/PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this UserObject"}>>> = 0
    phase<<<{"description": "The index corresponding to the phase for this UserObject"}>>> = 1
    flux_limiter_type<<<{"description": "Type of flux limiter to use.  'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = superbee
  []
  [afc_component1_phase0]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent<<<{"description": "Computes the advective flux of fluid of given phase and component.  Hence this UserObject is relevant to multi-phase, multi-component situations.  Explicitly, the UserObject computes (mass_fraction * density * relative_permeability / viscosity) * (- permeability * (grad(P) - density * gravity)), using the Kuzmin-Turek FEM-TVD multidimensional stabilization scheme", "href": "../../source/userobjects/PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this UserObject"}>>> = 1
    phase<<<{"description": "The index corresponding to the phase for this UserObject"}>>> = 0
    flux_limiter_type<<<{"description": "Type of flux limiter to use.  'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = superbee
  []
  [afc_component1_phase1]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent<<<{"description": "Computes the advective flux of fluid of given phase and component.  Hence this UserObject is relevant to multi-phase, multi-component situations.  Explicitly, the UserObject computes (mass_fraction * density * relative_permeability / viscosity) * (- permeability * (grad(P) - density * gravity)), using the Kuzmin-Turek FEM-TVD multidimensional stabilization scheme", "href": "../../source/userobjects/PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this UserObject"}>>> = 1
    phase<<<{"description": "The index corresponding to the phase for this UserObject"}>>> = 1
    flux_limiter_type<<<{"description": "Type of flux limiter to use.  'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = superbee
  []
[]

[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
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-3
  []
  [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)"}>>> = 2e7
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1e-5
  []
[]

[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
  []
  [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
  []
  [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-15 0 0 0 1e-15 0 0 0 1e-15'
  []
  [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"}>>>]
  [leftwater]
    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"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 3e6
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
  []
  [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"}>>> = 2e6
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = ppwater
  []
[]

[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_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'bcgs bjacobi 1E-15 1E-20 10000'
  []
[]

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

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

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  file_base<<<{"description": "Common file base name to be utilized with all output objects"}>>> = pressure_pulse_1d_2phasePS_KT
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false
  [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
  []
[]

and the UserObjects will look like

  [afc_component0_phase0]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 0
    phase = 0
    flux_limiter_type = superbee
  []
  [afc_component0_phase1]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 0
    phase = 1
    flux_limiter_type = superbee
  []
  [afc_component1_phase0]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 1
    phase = 0
    flux_limiter_type = superbee
  []
  [afc_component1_phase1]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 1
    phase = 1
    flux_limiter_type = superbee
  []
[]

[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 = 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
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 1
  []
[]

[BCs]
  [leftwater]
    type = DirichletBC
    boundary = left
    value = 3e6
    variable = ppwater
  []
  [rightwater]
    type = DirichletBC
    boundary = right
    value = 2e6
    variable = ppwater
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-15 1E-20 10000'
  []
[]

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

[VectorPostprocessors]
  [pp]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    sort_by = x
    variable = 'ppwater ppgas'
    start_point = '0 0 0'
    end_point = '100 0 0'
    num_points = 11
  []
[]

[Outputs]
  file_base = pressure_pulse_1d_2phasePS_KT
  print_linear_residuals = false
  [csv]
    type = CSV
    execute_on = final
  []
[]

References

D. Kuzmin and S. Turek. High-resolution FEM-TVD shcemes based on a fully multidimensional flux limiter. Journal of Computational Physics, 198:131–158, 2004.

@article{KuzminTurek2004,
    author = "Kuzmin, D. and Turek, S.",
    title = "{High-resolution FEM-TVD shcemes based on a fully multidimensional flux limiter}",
    journal = "Journal of Computational Physics",
    volume = "198",
    pages = "131--158",
    year = "2004"
}

(moose/modules/porous_flow/examples/tutorial/06_KT.i)

# Darcy flow with a tracer
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
  []
  [tracer_concentration]
  []
[]

[ICs]
  [tracer_concentration]
    type = FunctionIC
    function = '0.5*if(x*x+y*y<1.01,1,0)'
    variable = tracer_concentration
  []
[]

[PorousFlowFullySaturated]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  mass_fraction_vars = tracer_concentration
  stabilization = KT
  flux_limiter_type = superbee
[]

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1E6
    boundary = injection_area
  []
  [constant_outer_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 0
    boundary = rmax
  []
  [injected_tracer]
    type = DirichletBC
    variable = tracer_concentration
    value = 0.5
    boundary = injection_area
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2E9
    viscosity = 1.0E-3
    density0 = 1000.0
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_rel_tol = 1E-14
[]

[Outputs]
  exodus = true
[]

(moose/modules/porous_flow/examples/tutorial/08_KT.i)

# Unsaturated Darcy-Richards flow
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    input = annular
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
  []
[]

[PorousFlowUnsaturated]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  relative_permeability_exponent = 3
  relative_permeability_type = Corey
  residual_saturation = 0.1
  van_genuchten_alpha = 1E-6
  van_genuchten_m = 0.6
  stabilization = KT
  flux_limiter_type = None
[]

[BCs]
  [production]
    type = PorousFlowSink
    variable = porepressure
    fluid_phase = 0
    flux_function = 1E-2
    use_relperm = true
    boundary = injection_area
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2E9
    viscosity = 1.0E-3
    density0 = 1000.0
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E5
  dt = 1E5
  nl_rel_tol = 1E-14
[]

[Outputs]
  exodus = true
[]

(moose/modules/porous_flow/test/tests/heat_advection/heat_advection_1d.i)

# 1phase, heat advecting with a moving fluid
# Full upwinding is used
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 50
  xmin = 0
  xmax = 1
[]

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

[Variables]
  [temp]
    initial_condition = 200
  []
  [pp]
  []
[]

[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = '1-x'
  []
[]

[BCs]
  [pp0]
    type = DirichletBC
    variable = pp
    boundary = left
    value = 1
  []
  [pp1]
    type = DirichletBC
    variable = pp
    boundary = right
    value = 0
  []
  [spit_heat]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 300
  []
  [suck_heat]
    type = DirichletBC
    variable = temp
    boundary = right
    value = 200
  []
[]

[Kernels]
  [mass_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [advection]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pp
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = temp
  []
  [heat_advection]
    type = PorousFlowHeatAdvection
    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 = PorousFlowTemperature
    temperature = temp
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1.0
    density = 125
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 2 0 0 0 3'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [PS]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres bjacobi 1E-15 1E-10 10000'
  []
[]

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

[VectorPostprocessors]
  [T]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '1 0 0'
    num_points = 51
    sort_by = x
    variable = temp
  []
[]

[Outputs]
  [csv]
    type = CSV
    sync_times = '0.1 0.6'
    sync_only = true
  []
[]

(moose/modules/porous_flow/test/tests/heat_advection/heat_advection_1d_KT.i)

# 1phase, heat advecting with a moving fluid
# Using the Kuzmin-Turek stabilization scheme
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 50
  xmin = 0
  xmax = 1
[]

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

[Variables]
  [temp]
    initial_condition = 200
  []
  [pp]
  []
[]

[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = '1-x'
  []
[]

[BCs]
  [pp0]
    type = DirichletBC
    variable = pp
    boundary = left
    value = 1
  []
  [pp1]
    type = DirichletBC
    variable = pp
    boundary = right
    value = 0
  []
  [spit_heat]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 300
  []
  [suck_heat]
    type = DirichletBC
    variable = temp
    boundary = right
    value = 200
  []
[]

[Kernels]
  [mass_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [fluid_advection]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = pp
    advective_flux_calculator = fluid_advective_flux
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = temp
  []
  [heat_advection]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = temp
    advective_flux_calculator = heat_advective_flux
  []
[]

[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
  []
  [fluid_advective_flux]
    type = PorousFlowAdvectiveFluxCalculatorSaturated
    flux_limiter_type = superbee
  []
  [heat_advective_flux]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedHeat
    flux_limiter_type = superbee
  []
[]

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

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1.0
    density = 125
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 2 0 0 0 3'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [PS]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres bjacobi 1E-15 1E-10 10000'
  []
[]

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

[VectorPostprocessors]
  [T]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '1 0 0'
    num_points = 51
    sort_by = x
    variable = temp
  []
[]

[Outputs]
  file_base = heat_advection_1d_KT
  [csv]
    type = CSV
    sync_times = '0.1 0.6'
    sync_only = true
  []
[]

(moose/modules/porous_flow/test/tests/heat_advection/heat_advection_1d_KT.i)

# 1phase, heat advecting with a moving fluid
# Using the Kuzmin-Turek stabilization scheme
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 50
  xmin = 0
  xmax = 1
[]

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

[Variables]
  [temp]
    initial_condition = 200
  []
  [pp]
  []
[]

[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = '1-x'
  []
[]

[BCs]
  [pp0]
    type = DirichletBC
    variable = pp
    boundary = left
    value = 1
  []
  [pp1]
    type = DirichletBC
    variable = pp
    boundary = right
    value = 0
  []
  [spit_heat]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 300
  []
  [suck_heat]
    type = DirichletBC
    variable = temp
    boundary = right
    value = 200
  []
[]

[Kernels]
  [mass_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [fluid_advection]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = pp
    advective_flux_calculator = fluid_advective_flux
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = temp
  []
  [heat_advection]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = temp
    advective_flux_calculator = heat_advective_flux
  []
[]

[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
  []
  [fluid_advective_flux]
    type = PorousFlowAdvectiveFluxCalculatorSaturated
    flux_limiter_type = superbee
  []
  [heat_advective_flux]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedHeat
    flux_limiter_type = superbee
  []
[]

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

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1.0
    density = 125
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 2 0 0 0 3'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [PS]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'gmres bjacobi 1E-15 1E-10 10000'
  []
[]

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

[VectorPostprocessors]
  [T]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '1 0 0'
    num_points = 51
    sort_by = x
    variable = temp
  []
[]

[Outputs]
  file_base = heat_advection_1d_KT
  [csv]
    type = CSV
    sync_times = '0.1 0.6'
    sync_only = true
  []
[]

(moose/modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D.i)

# Using flux-limited TVD advection ala Kuzmin and Turek, mploying PorousFlow Kernels and UserObjects, with superbee flux-limiter
# 1D version
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmin = 0
  xmax = 1
[]

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

[Variables]
  [porepressure]
  []
  [tracer]
  []
[]

[ICs]
  [porepressure]
    type = FunctionIC
    variable = porepressure
    function = '1 - x'
  []
  [tracer]
    type = FunctionIC
    variable = tracer
    function = 'if(x<0.1,0,if(x>0.3,0,1))'
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = tracer
  []
  [flux0]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = tracer
    advective_flux_calculator = advective_flux_calculator_0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = porepressure
  []
  [flux1]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = porepressure
    advective_flux_calculator = advective_flux_calculator_1
  []
[]

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1
    boundary = left
  []
  [no_tracer_on_left]
    type = DirichletBC
    variable = tracer
    value = 0
    boundary = left
  []
  [remove_component_1]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    boundary = right
    fluid_phase = 0
    pt_vals = '0 1E3'
    multipliers = '0 1E3'
    mass_fraction_component = 1
    use_mobility = true
    flux_function = 1E3
  []
  [remove_component_0]
    type = PorousFlowPiecewiseLinearSink
    variable = tracer
    boundary = right
    fluid_phase = 0
    pt_vals = '0 1E3'
    multipliers = '0 1E3'
    mass_fraction_component = 0
    use_mobility = true
    flux_function = 1E3
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2E9
    thermal_expansion = 0
    viscosity = 1.0
    density0 = 1000.0
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure tracer'
    number_fluid_phases = 1
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
  []
  [advective_flux_calculator_0]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
    flux_limiter_type = superbee
    fluid_component = 0
  []
  [advective_flux_calculator_1]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
    flux_limiter_type = superbee
    fluid_component = 1
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = porepressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = tracer
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = the_simple_fluid
    phase = 0
  []
  [relperm]
    type = PorousFlowRelativePermeabilityConst
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-2 0 0   0 1E-2 0   0 0 1E-2'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[VectorPostprocessors]
  [tracer]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '1 0 0'
    num_points = 11
    sort_by = x
    variable = tracer
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 6
  dt = 6E-2
  nl_abs_tol = 1E-8
  timestep_tolerance = 1E-3
[]

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

(moose/modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D.i)

# Using flux-limited TVD advection ala Kuzmin and Turek, mploying PorousFlow Kernels and UserObjects, with superbee flux-limiter
# 1D version
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmin = 0
  xmax = 1
[]

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

[Variables]
  [porepressure]
  []
  [tracer]
  []
[]

[ICs]
  [porepressure]
    type = FunctionIC
    variable = porepressure
    function = '1 - x'
  []
  [tracer]
    type = FunctionIC
    variable = tracer
    function = 'if(x<0.1,0,if(x>0.3,0,1))'
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = tracer
  []
  [flux0]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = tracer
    advective_flux_calculator = advective_flux_calculator_0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = porepressure
  []
  [flux1]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = porepressure
    advective_flux_calculator = advective_flux_calculator_1
  []
[]

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1
    boundary = left
  []
  [no_tracer_on_left]
    type = DirichletBC
    variable = tracer
    value = 0
    boundary = left
  []
  [remove_component_1]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    boundary = right
    fluid_phase = 0
    pt_vals = '0 1E3'
    multipliers = '0 1E3'
    mass_fraction_component = 1
    use_mobility = true
    flux_function = 1E3
  []
  [remove_component_0]
    type = PorousFlowPiecewiseLinearSink
    variable = tracer
    boundary = right
    fluid_phase = 0
    pt_vals = '0 1E3'
    multipliers = '0 1E3'
    mass_fraction_component = 0
    use_mobility = true
    flux_function = 1E3
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2E9
    thermal_expansion = 0
    viscosity = 1.0
    density0 = 1000.0
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure tracer'
    number_fluid_phases = 1
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
  []
  [advective_flux_calculator_0]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
    flux_limiter_type = superbee
    fluid_component = 0
  []
  [advective_flux_calculator_1]
    type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
    flux_limiter_type = superbee
    fluid_component = 1
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = porepressure
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = tracer
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = the_simple_fluid
    phase = 0
  []
  [relperm]
    type = PorousFlowRelativePermeabilityConst
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-2 0 0   0 1E-2 0   0 0 1E-2'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[VectorPostprocessors]
  [tracer]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '1 0 0'
    num_points = 11
    sort_by = x
    variable = tracer
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 6
  dt = 6E-2
  nl_abs_tol = 1E-8
  timestep_tolerance = 1E-3
[]

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

(moose/modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_2phasePS_KT.i)

# Pressure pulse in 1D with 2 phases, 2components - transient
# Using KT stabilization

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

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

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

[AuxVariables]
  [massfrac_ph0_sp0]
    initial_condition = 1
  []
  [massfrac_ph1_sp0]
    initial_condition = 0
  []
  [ppgas]
    family = MONOMIAL
    order = FIRST
  []
[]

[Kernels]
  [mass_component0]
    type = PorousFlowMassTimeDerivative
    variable = ppwater
    fluid_component = 0
  []
  [flux_component0_phase0]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = ppwater
    advective_flux_calculator = afc_component0_phase0
  []
  [flux_component0_phase1]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = ppwater
    advective_flux_calculator = afc_component0_phase1
  []
  [mass_component1]
    type = PorousFlowMassTimeDerivative
    variable = sgas
    fluid_component = 1
  []
  [flux_component1_phase0]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = sgas
    advective_flux_calculator = afc_component1_phase0
  []
  [flux_component1_phase1]
    type = PorousFlowFluxLimitedTVDAdvection
    variable = sgas
    advective_flux_calculator = afc_component1_phase1
  []
[]

[AuxKernels]
  [ppgas]
    type = PorousFlowPropertyAux
    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
  []
  [afc_component0_phase0]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 0
    phase = 0
    flux_limiter_type = superbee
  []
  [afc_component0_phase1]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 0
    phase = 1
    flux_limiter_type = superbee
  []
  [afc_component1_phase0]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 1
    phase = 0
    flux_limiter_type = superbee
  []
  [afc_component1_phase1]
    type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
    fluid_component = 1
    phase = 1
    flux_limiter_type = superbee
  []
[]

[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 = 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
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 0
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 1
    phase = 1
  []
[]

[BCs]
  [leftwater]
    type = DirichletBC
    boundary = left
    value = 3e6
    variable = ppwater
  []
  [rightwater]
    type = DirichletBC
    boundary = right
    value = 2e6
    variable = ppwater
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value = 'bcgs bjacobi 1E-15 1E-20 10000'
  []
[]

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

[VectorPostprocessors]
  [pp]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    sort_by = x
    variable = 'ppwater ppgas'
    start_point = '0 0 0'
    end_point = '100 0 0'
    num_points = 11
  []
[]

[Outputs]
  file_base = pressure_pulse_1d_2phasePS_KT
  print_linear_residuals = false
  [csv]
    type = CSV
    execute_on = final
  []
[]