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.

commentnote

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]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  mass_fraction_vars = tracer_concentration
  stabilization = KT
  flux_limiter_type = superbee
[]
(modules/porous_flow/examples/tutorial/06_KT.i)

and

[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
[]
(modules/porous_flow/examples/tutorial/08_KT.i)

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]
  [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
  []
[]
(modules/porous_flow/test/tests/heat_advection/heat_advection_1d.i)

while an identical simulation using KT stabilization has Kernels:

[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
  []
[]
(modules/porous_flow/test/tests/heat_advection/heat_advection_1d_KT.i)

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

EquationUserObjectUse case
PorousFlowAdvectiveFluxCalculatorSaturated1-phase, 1-component, fully-saturated
PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent1-phase, multi-component, fully-saturated
PorousFlowAdvectiveFluxCalculatorUnsaturated1-phase, 1-component, unsaturated
PorousFlowAdvectiveFluxCalculatorUnsaturatedmulti-phase, multi-component, immiscible
PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponentmulti-phase, multi-component
PorousFlowAdvectiveFluxCalculatorSaturatedHeatheat, 1-phase
PorousFlowAdvectiveFluxCalculatorUnsaturatedHeatheat, multi-phase

The notation is detailed in the nomenclature: briefly, a subscript denotes a phase, a superscript denotes a fluid component, is a fluid density, a fluid viscosity, the permeability tensor, a fluid porepressure, the acceleration of gravity, a mass fraction, a relative permeability, and 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
  []
[]
(modules/porous_flow/test/tests/heat_advection/heat_advection_1d_KT.i)

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]
  [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
  []
[]
(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D.i)
  [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
  []
[]
(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D.i)
commentnote

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]
  [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
  []
[]
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_2phasePS_KT.i)

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
  []
[]
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_2phasePS_KT.i)

References

  1. 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.[BibTeX]