# 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
[../]
fluid_component = 0
variable = pp
[../]
[./energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temp
[../]
variable = temp
[../]


while an identical simulation using KT stabilization has Kernels:

[Kernels]
[./mass_dot]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[../]
variable = pp
[../]
[./energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temp
[../]
variable = temp
[../]


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

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]
flux_limiter_type = superbee
[../]
flux_limiter_type = superbee
[../]


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]
variable = tracer
[../]
[./mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[../]
[./flux1]
variable = porepressure
[../]
[]

(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D.i)
  [./advective_flux_calculator_0]
flux_limiter_type = superbee
fluid_component = 0
[../]
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]
variable = ppwater
[../]
[./flux_component0_phase1]
variable = ppwater
[../]
[./mass_component1]
type = PorousFlowMassTimeDerivative
variable = sgas
fluid_component = 1
[../]
[./flux_component1_phase0]
variable = sgas
[../]
[./flux_component1_phase1]
variable = sgas
[../]
[]

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

and the UserObjects will look like

  [./afc_component0_phase0]
fluid_component = 0
phase = 0
flux_limiter_type = superbee
[../]
[./afc_component0_phase1]
fluid_component = 0
phase = 1
flux_limiter_type = superbee
[../]
[./afc_component1_phase0]