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