Numerical diffusion

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

Numerical diffusion is the artificial smoothing of quantities, such as temperature and concentrations, as they are transported through a numerical model. In the case of PorousFlow, it is usually the fluid that transports these quantities (the fluid "advects" temperature and chemical species). Numerical diffusion can be the major source of inaccurate results in simulations, as MOOSE predicts that tracer breakthrough times (etc) are much shorter than they are in reality.

An animation of the tracer advection is shown in Figure 1. Notice that the initially-sharp profile suffers from diffusion.

Figure 1: Tracer advection down the porepressure gradient.

Numerical diffusion has been mentioned in various pieces of PorousFlow documentation (eg, in the Porous Flow Tutorial Page 06. Adding a tracer and PorousFlowFullySaturatedDarcyFlow, and the documentation concerning the sinks tests and heat advection test page). There are two main sources of numerical diffusion:

  • employing large time steps with MOOSE's implicit time stepping scheme

  • the full upwinding used by PorousFlow

To quantify the numerical diffusion by example, a single-phase tracer advection problem is studied in 1D. A porepressure gradient is established so that the Darcy velocity, m/s, and porosity is chosen to be so that the tracer advects down the porepressure gradient with a constant velocity of m/s.

The degree of numerical diffusion is a function of the spatial and temporal discretisation, as well as the upwinding, RDG reconstruction and limiting. To quantify this, seven input files are created:

  1. Using framework MOOSE objects: there is no mass lumping and no upwinding.

  2. Using the PorousFlow fully-saturated action, which is mass-lumped but does not use any upwinding.

  3. Using the standard mass-lumped and fully-upwinded PorousFlow kernels.

  4. Using the RDG module with no reconstruction (RDG(P0)): the tracer is a constant monomial Variable (c.f. linear lagrange variables in PorousFlow)

  5. Using the RDG module with linear reconstruction (RDG(P0P1)): the tracer is constant monomial, but a linear form is reconstructed and then limited using RDG's mc limiter . Another limiter could have been chosen, with very similar results.

  6. Using the Kuzmin-Turek scheme with no flux limiter (should be identical to the fully upwinded case).

  7. Using the Kuzmin-Turek scheme with a flux limiter (should be similar to RDG(POP1) case).

The Kuzmin-Turek scheme is described in: D Kuzmin and S Turek (2004) "High-resolution FEM-TVD schemes based on a fully multidimensional flux limiter." Journal of Computational Physics, Volume 198, 131 - 158.

No mass lumping and no upwinding

The MOOSE input file is:

# Using framework objects: no mass lumping or upwinding
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmin = 0
  xmax = 1
[]

[Variables]
  [tracer]
  []
[]

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

[Kernels]
  [mass_dot]
    type = TimeDerivative
    variable = tracer
  []
  [flux]
    type = ConservativeAdvection
    velocity = '0.1 0 0'
    variable = tracer
  []
[]

[BCs]
  [no_tracer_on_left]
    type = DirichletBC
    variable = tracer
    value = 0
    boundary = left
  []
  [remove_tracer]
    # Ideally, an OutflowBC would be used, but that does not exist in the framework
    # In 1D VacuumBC is the same as OutflowBC, with the alpha parameter being twice the velocity
    type = VacuumBC
    boundary = right
    alpha = 0.2 # 2 * velocity
    variable = tracer
  []
[]

[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 = 101
    sort_by = x
    variable = tracer
  []
[]

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

[Outputs]
  [out]
    type = CSV
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/numerical_diffusion/framework.i)

Figure 2 and Figure 3 show the dependence on discretisation when there is no mass-lumping and no upwinding. Evidently, the lack of upwinding causes overshoots and undershoots, and the diffusion becomes less as the spatial discretisation becomes finer.

Figure 2: Diffusion as a function of number of elements. The number of time steps is fixed to 100 and no mass lumping and no upwinding is used.

Figure 3: Diffusion as a function of number of time steps. The number of elements is fixed to 100 and no mass lumping and no upwinding is used.

Mass lumping and no upwinding

The MOOSE input file is:

# Using the fully-saturated action, which does mass lumping but no upwinding
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmin = 0
  xmax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[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))'
  []
[]

[PorousFlowFullySaturated]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  mass_fraction_vars = tracer
  stabilization = none
[]

[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
  []
[]

[Materials]
  [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 = 101
    sort_by = x
    variable = tracer
  []
[]

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

[Outputs]
  [out]
    type = CSV
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/numerical_diffusion/fully_saturated_action.i)

Figure 4 and Figure 5 show the dependence on discretisation when there is no upwinding. Evidently, the lack of upwinding causes overshoots and undershoots, and the diffusion becomes less as the spatial discretisation becomes finer.

Figure 4: Diffusion as a function of number of elements. The number of time steps is fixed to 100 and no upwinding is used.

Figure 5: Diffusion as a function of number of time steps. The number of elements is fixed to 100 and no upwinding is used.

Mass lumping and full upwinding

The MOOSE input file is identical to the one above, save for the PorousFlowFullySaturated block that must now contain stabilization = Full:


[PorousFlowFullySaturated]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  mass_fraction_vars = tracer
  stabilization = Full
[]

(Also, the test suite contains a non-action version).

Figure 6 and Figure 7 show the dependence on discretisation when full upwinding is used.

Figure 6: Diffusion as a function of number of elements. The number of time steps is fixed to 100 and full upwinding is used.

Figure 7: Diffusion as a function of number of time steps. The number of elements is fixed to 100 and full upwinding is used.

RDG(P0)

In RDG, the variables (just the tracer in this case) are constant-monomial, that is, they are constant throughout the element. This means MOOSE behaves like a Finite Volume Method, and it means that the output looks less smooth than the usual linear-Lagrange variables (the figures below display "stepped" results). The MOOSE input file is:

# This test demonstrates the advection of a tracer in 1D using the RDG module.
# There is no slope limiting.  Changing the SlopeLimiting scheme to minmod, mc,
# or superbee means that a linear reconstruction is performed, and the slope
# limited according to the scheme chosen.  Doing this produces RDG(P0P1) and
# substantially reduces numerical diffusion

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

[Variables]
  [./tracer]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

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

[UserObjects]
  [./lslope]
    type = AEFVSlopeLimitingOneD
    execute_on = 'linear'
    scheme = 'none' #none | minmod | mc | superbee
    u = tracer
  [../]

  [./internal_side_flux]
    type = AEFVUpwindInternalSideFlux
    execute_on = 'linear'
    velocity = 0.1
  [../]

  [./free_outflow_bc]
    type = AEFVFreeOutflowBoundaryFlux
    execute_on = 'linear'
    velocity = 0.1
  [../]
[]

[Kernels]
  [./dot]
    type = TimeDerivative
    variable = tracer
  [../]
[]
[DGKernels]
  [./concentration]
    type = AEFVKernel
    variable = tracer
    component = 'concentration'
    flux = internal_side_flux
    u = tracer
  [../]
[]

[BCs]
  [./concentration]
    type = AEFVBC
    boundary = 'left right'
    variable = tracer
    component = 'concentration'
    flux = free_outflow_bc
    u = tracer
  [../]
[]

[Materials]
  [./aefv]
    type = AEFVMaterial
    slope_limiting = lslope
    u = tracer
  [../]
[]

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

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

[Outputs]
  #exodus = true
  csv = true
  execute_on = final
[]
(modules/rdg/test/tests/advection_1d/rdgP0.i)

Figure 8 and Figure 9 show the dependence on discretisation when RDG(P0) is used. This is identical to mass-lumping + full-upwinding (up to the fact that constant-monomial variables are used instead of linear-Lagrange). As expected, there are no oscillations or over-shoots or under-shoots, but the results suffer from numerical diffusion.

Figure 8: Diffusion as a function of number of elements. The number of time steps is fixed to 100 and RDG(P0) is used.

Figure 9: Diffusion as a function of number of time steps. The number of elements is fixed to 100 and RDG(P0) is used.

RDG(P0P1)

The MOOSE input file is almost identical to the RDG(P0) version, save for the changing the SlopeLimiting scheme from "none" to "mc" (or another limiter such as "minmod" or "superbee").

Figure 10 and Figure 11 show the dependence on discretisation when RDG(P0P1) is used. As with RDG(P0), there are no oscillations or over-shoots or under-shoots. However, numerical diffusion is greatly reduced by the linear reconstruction.

Figure 10: Diffusion as a function of number of elements. The number of time steps is fixed to 100 and RDG(P0P1) is used.

Figure 11: Diffusion as a function of number of time steps. The number of elements is fixed to 100 and RDG(P0P1) is used.

Kuzmin-Turek scheme with no flux limiter

To employ the Kuzmin-Turek scheme, the fully_saturated_action.i input file listed above may be simply modified by including KT stabilization:


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

in the PorousFlowFullySaturated Action block. Alternatively, a new input file may be built, such as:

# Using Flux-Limited TVD Advection ala Kuzmin and Turek, but without any antidiffusion
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmin = 0
  xmax = 1
[]

[Variables]
  [tracer]
  []
[]

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

[Kernels]
  [mass_dot]
    type = MassLumpedTimeDerivative
    variable = tracer
  []
  [flux]
    type = FluxLimitedTVDAdvection
    variable = tracer
    advective_flux_calculator = fluo
  []
[]

[UserObjects]
  [fluo]
    type = AdvectiveFluxCalculatorConstantVelocity
    flux_limiter_type = none
    u = tracer
    velocity = '0.1 0 0'
  []
[]

[BCs]
  [no_tracer_on_left]
    type = DirichletBC
    variable = tracer
    value = 0
    boundary = left
  []
  [remove_tracer]
# Ideally, an OutflowBC would be used, but that does not exist in the framework
# In 1D VacuumBC is the same as OutflowBC, with the alpha parameter being twice the velocity
    type = VacuumBC
    boundary = right
    alpha = 0.2 # 2 * velocity
    variable = tracer
  []
[]

[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 = 101
    sort_by = x
    variable = tracer
  []
[]

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

[Outputs]
  [out]
    type = CSV
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/numerical_diffusion/fltvd_no_antidiffusion.i)

Figure 12 and Figure 13 show the dependence on discretisation when Kuzmin-Turek with no flux limiter is used. As expected, the behaviour is identical to the case with full upwinding.

Figure 12: Diffusion as a function of number of elements. The number of time steps is fixed to 100 and Kuzmin-Turek with no flux limiter is used.

Figure 13: Diffusion as a function of number of time steps. The number of elements is fixed to 100 and Kuzmin-Turek with no flux limiter is used.

Kuzmin-Turek scheme with flux limiter

To employ the Kuzmin-Turek scheme, the fully_saturated_action.i input file listed above may be simply modified by including KT stabilization:


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

in the PorousFlowFullySaturated Action block. Alternatively, a new input file may be built, such as:

# Using Flux-Limited TVD Advection ala Kuzmin and Turek, with antidiffusion from superbee flux limiting
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmin = 0
  xmax = 1
[]

[Variables]
  [tracer]
  []
[]

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

[Kernels]
  [mass_dot]
    type = MassLumpedTimeDerivative
    variable = tracer
  []
  [flux]
    type = FluxLimitedTVDAdvection
    variable = tracer
    advective_flux_calculator = fluo
  []
[]

[UserObjects]
  [fluo]
    type = AdvectiveFluxCalculatorConstantVelocity
    flux_limiter_type = superbee
    u = tracer
    velocity = '0.1 0 0'
  []
[]

[BCs]
  [no_tracer_on_left]
    type = DirichletBC
    variable = tracer
    value = 0
    boundary = left
  []
  [remove_tracer]
# Ideally, an OutflowBC would be used, but that does not exist in the framework
# In 1D VacuumBC is the same as OutflowBC, with the alpha parameter being twice the velocity
    type = VacuumBC
    boundary = right
    alpha = 0.2 # 2 * velocity
    variable = tracer
  []
[]

[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 = 101
    sort_by = x
    variable = tracer
  []
[]

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

[Outputs]
  [out]
    type = CSV
    execute_on = final
  []
[]
(modules/porous_flow/test/tests/numerical_diffusion/fltvd.i)

Figure 14 and Figure 15 show the dependence on discretisation when Kuzmin-Turek with a flux limiter is used (the "superbee" limiter has been chosen in this case, see the worked example of Kuzmin-Turek stabilization). As expected, the behaviour is similar to the RDG(POP1) case, except it is smooth because the variables are linear-lagrange.

Figure 14: Diffusion as a function of number of elements. The number of time steps is fixed to 100 and Kuzmin-Turek with the superbee flux limiter.

Figure 15: Diffusion as a function of number of time steps. The number of elements is fixed to 100 and Kuzmin-Turek with the superbee flux limiter.

Summary

By way of comparison, Figure 16 and Figure 17 show the results when the number of elements is fixed to 100 and the number of timesteps is 100 or 1000 (respectively). Evidently, without upwinding, the result includes over-shoots and under-shoots, which are typically devastating for PorousFlow simulations, since they typically manifest themselves as negative saturations, negative mass fractions, or negative temperatures. With upwinding, RDG(P0), or Kuzmin-Turek with no antidiffusion these are avoided, at the expense of introducing large numerical diffusion. RDG(P0P1) and Kuzmin-Turek with a flux limiter produces a similar amount of numerical diffusion to the non-upwinded, non-masslumped version, but without introducing any over-shoots and under-shoots.

Figure 16: Diffusion for different numerical schemes. 100 timesteps are used.

Figure 17: Diffusion for different numerical schemes. 1000 timesteps are used.

Interestingly, as the timestep size is decreased, the numerical diffusion is decreased when using lumping and upwinding, but only up to a point: increasing the number of timesteps from 100 to 10000 makes basically no difference when the number of elements is 100. This suggests that lumping and full-upwinding will always produce numerical diffusion, and this may be proved by simple analysis (not presented here) since the tracer always gets moved from upwind node to downwind node irrespective of the concentration of the tracer.