Sinks test descriptions

Introduction

As described in the page on boundary conditions, a number of different sink boundary conditions have been implemented in PorousFlow. This page describes some of the tests of the various forms of sink/source boundary conditions.

To make these into sources instead of sinks, the strength of the flux just needs to be made negative. All the sinks are implemented using full upwinding. This is to prevent the sink from attempting to remove fluid from a node that actually contains no fluid.

The basic sink uses a Function to specify the flux on the boundary, and also has the option of multiplying by any combination of: the fluid mobility, the relative permeability, or a mass fraction. These latter multiplying factors are all useful in the case of sinks to prevent an unlimited amount of fluid being withdrawn from the porous medium, which can lead to extremely poor nonlinear convergence even if only one node in the entire mesh is "running dry".

Derived from the basic one, is another boundary condition that allows the flux to be modified by a piecewise-linear function of porepressure, which is useful for the case where transfer coefficients are defined across the boundary, or more complicated situations.

Also derived from the basic one are two others, in which the flux is governed by a half Gaussian or half cubic function of porepressure, which are useful for modelling evapotranspiration through a boundary.

Basic PorousFlow Sink

Test 1

A sink flux of strength 6kg.m.s is applied to the left edge () of a 3D mesh. A single-phase, single-component fluid is used, and the porepressure is initialised to (for ). No fluid flow within the element is used, so the masses of fluid at the finite-element nodes behave independently. The fluid is assumed to have density kg.m. The porosity is 0.1.

The input file:

# apply a sink flux and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  [../]
[]

[Variables]
  [./pp]
  [../]
[]

[ICs]
  [./pp]
    type = FunctionIC
    variable = pp
    function = y+1
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-5 0 0 0 1E-5 0 0 0 1E-5'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 0
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
  [./xval]
  [../]
  [./yval]
  [../]
[]

[ICs]
  [./xval]
    type = FunctionIC
    variable = xval
    function = x
  [../]
  [./yval]
    type = FunctionIC
    variable = yval
    function = y
  [../]
[]

[Functions]
  [./mass00]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p00 1.3'
  [../]
  [./mass01]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p01 1.3'
  [../]
  [./expected_mass_change00]
    type = ParsedFunction
    value = 'fcn*perm*dens0*exp(pp/bulk)/visc*area*dt'
    vars = 'fcn perm dens0 pp bulk visc area dt'
    vals = '6   1    1      0  1.3  1  0.5  1E-3'
  [../]
[]

[Postprocessors]
  [./p00]
    type = PointValue
    point = '0 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m00]
    type = FunctionValuePostprocessor
    function = mass00
    execute_on = 'initial timestep_end'
  [../]
  [./del_m00]
    type = FunctionValuePostprocessor
    function = expected_mass_change00
    execute_on = 'timestep_end'
  [../]
  [./p10]
    type = PointValue
    point = '1 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p01]
    type = PointValue
    point = '0 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m01]
    type = FunctionValuePostprocessor
    function = mass01
    execute_on = 'initial timestep_end'
  [../]
  [./p11]
    type = PointValue
    point = '1 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
[]

[BCs]
  [./flux]
    type = PorousFlowSink
    boundary = 'left'
    variable = pp
    use_mobility = false
    use_relperm = true
    fluid_phase = 0
    flux_function = 6
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-3
  end_time = 1E-2
  nl_rel_tol = 1E-12
  nl_abs_tol = 1E-12
[]

[Outputs]
  file_base = s01
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
  [../]
  [./csv]
    type = CSV
    execute_on = 'initial timestep_end'
  [../]
[]
(modules/porous_flow/test/tests/sinks/s01.i)

Under these conditions, assuming so that the porous medium is fully saturated with fluid, the fluid mass at a node should obey (1) where is the volume occupied by the node, and is its area exposed to the flux. MOOSE correctly produces this result, as illustrated in Figure 1.

Figure 1: Results of Test 1, illustrating that MOOSE correctly applies a constant sink flux to boundary nodes.

Test 2

An identical setup to Test 1 is used here, but with the sink flux strength being multiplied by the mobility: (2) where is the permeability tensor projected onto the normal direction to the boundary , and the fluid density and viscosity are and , respectively. In this example Pa.s and m. The other parameters are the same as Test 1, except now the strength of the flux is 6Pa.s.

The input file:

# apply a sink flux with use_mobility=true and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  [../]
[]

[Variables]
  [./pp]
  [../]
[]

[ICs]
  [./pp]
    type = FunctionIC
    variable = pp
    function = y+1
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.2 0 0 0 0.1 0 0 0 0.1'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 0
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
  [./xval]
  [../]
  [./yval]
  [../]
[]

[ICs]
  [./xval]
    type = FunctionIC
    variable = xval
    function = x
  [../]
  [./yval]
    type = FunctionIC
    variable = yval
    function = y
  [../]
[]

[Functions]
  [./mass00]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p00 1.3'
  [../]
  [./mass01]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p01 1.3'
  [../]
  [./expected_mass_change00]
    type = ParsedFunction
    value = 'fcn*perm*dens0*exp(pp/bulk)/visc*area*dt'
    vars = 'fcn perm dens0 pp bulk visc area dt'
    vals = '6   0.2  1.1  p00 1.3  1.1  0.5  1E-3'
  [../]
  [./expected_mass_change01]
    type = ParsedFunction
    value = 'fcn*perm*dens0*exp(pp/bulk)/visc*area*dt'
    vars = 'fcn perm dens0 pp bulk visc area dt'
    vals = '6   0.2  1.1  p01 1.3  1.1  0.5  1E-3'
  [../]
  [./mass00_expect]
    type = ParsedFunction
    value = 'mass_prev-mass_change'
    vars = 'mass_prev mass_change'
    vals = 'm00_prev  del_m00'
  [../]
  [./mass01_expect]
    type = ParsedFunction
    value = 'mass_prev-mass_change'
    vars = 'mass_prev mass_change'
    vals = 'm01_prev  del_m01'
  [../]
[]

[Postprocessors]
  [./p00]
    type = PointValue
    point = '0 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m00]
    type = FunctionValuePostprocessor
    function = mass00
    execute_on = 'initial timestep_end'
  [../]
  [./m00_prev]
    type = FunctionValuePostprocessor
    function = mass00
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./del_m00]
    type = FunctionValuePostprocessor
    function = expected_mass_change00
    execute_on = 'timestep_end'
    outputs = 'console'
  [../]
  [./m00_expect]
    type = FunctionValuePostprocessor
    function = mass00_expect
    execute_on = 'timestep_end'
  [../]
  [./p10]
    type = PointValue
    point = '1 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p01]
    type = PointValue
    point = '0 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m01]
    type = FunctionValuePostprocessor
    function = mass01
    execute_on = 'initial timestep_end'
  [../]
  [./m01_prev]
    type = FunctionValuePostprocessor
    function = mass01
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./del_m01]
    type = FunctionValuePostprocessor
    function = expected_mass_change01
    execute_on = 'timestep_end'
    outputs = 'console'
  [../]
  [./m01_expect]
    type = FunctionValuePostprocessor
    function = mass01_expect
    execute_on = 'timestep_end'
  [../]
  [./p11]
    type = PointValue
    point = '1 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
[]

[BCs]
  [./flux]
    type = PorousFlowSink
    boundary = 'left'
    variable = pp
    use_mobility = true
    use_relperm = true
    fluid_phase = 0
    flux_function = 6
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-3
  end_time = 0.03
  nl_rel_tol = 1E-12
  nl_abs_tol = 1E-12
[]

[Outputs]
  file_base = s02
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
    interval = 30
  [../]
  [./csv]
    type = CSV
    execute_on = 'timestep_end'
    interval = 3
  [../]
[]
(modules/porous_flow/test/tests/sinks/s02.i)

In this case, the expected result is (for ) (3) MOOSE correctly produces this result, as illustrated in Figure 2.

Figure 2: Results of Test 2, illustrating that MOOSE correctly applies a constant sink flux modified by the fluid mobility. (A slight drift away from the expected result is due to MOOSE taking large time steps.)

Test 3

An identical setup to Test 1 is used here, but with the sink flux strength being multiplied by the relative permeability, which is chosen to be: (4) with being the fluid saturation. A van-Genuchten capillary relationship is used: (5) with Pa, and . The porepressure is initialised to be . The other parameters are identical to Test 1.

The input file:

# apply a sink flux with use_mobility=true and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  [../]
[]

[Variables]
  [./pp]
  [../]
[]

[ICs]
  [./pp]
    type = FunctionIC
    variable = pp
    function = y+1
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.2 0 0 0 0.1 0 0 0 0.1'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 0
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
  [./xval]
  [../]
  [./yval]
  [../]
[]

[ICs]
  [./xval]
    type = FunctionIC
    variable = xval
    function = x
  [../]
  [./yval]
    type = FunctionIC
    variable = yval
    function = y
  [../]
[]

[Functions]
  [./mass00]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p00 1.3'
  [../]
  [./mass01]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p01 1.3'
  [../]
  [./expected_mass_change00]
    type = ParsedFunction
    value = 'fcn*perm*dens0*exp(pp/bulk)/visc*area*dt'
    vars = 'fcn perm dens0 pp bulk visc area dt'
    vals = '6   0.2  1.1  p00 1.3  1.1  0.5  1E-3'
  [../]
  [./expected_mass_change01]
    type = ParsedFunction
    value = 'fcn*perm*dens0*exp(pp/bulk)/visc*area*dt'
    vars = 'fcn perm dens0 pp bulk visc area dt'
    vals = '6   0.2  1.1  p01 1.3  1.1  0.5  1E-3'
  [../]
  [./mass00_expect]
    type = ParsedFunction
    value = 'mass_prev-mass_change'
    vars = 'mass_prev mass_change'
    vals = 'm00_prev  del_m00'
  [../]
  [./mass01_expect]
    type = ParsedFunction
    value = 'mass_prev-mass_change'
    vars = 'mass_prev mass_change'
    vals = 'm01_prev  del_m01'
  [../]
[]

[Postprocessors]
  [./p00]
    type = PointValue
    point = '0 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m00]
    type = FunctionValuePostprocessor
    function = mass00
    execute_on = 'initial timestep_end'
  [../]
  [./m00_prev]
    type = FunctionValuePostprocessor
    function = mass00
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./del_m00]
    type = FunctionValuePostprocessor
    function = expected_mass_change00
    execute_on = 'timestep_end'
    outputs = 'console'
  [../]
  [./m00_expect]
    type = FunctionValuePostprocessor
    function = mass00_expect
    execute_on = 'timestep_end'
  [../]
  [./p10]
    type = PointValue
    point = '1 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p01]
    type = PointValue
    point = '0 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m01]
    type = FunctionValuePostprocessor
    function = mass01
    execute_on = 'initial timestep_end'
  [../]
  [./m01_prev]
    type = FunctionValuePostprocessor
    function = mass01
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./del_m01]
    type = FunctionValuePostprocessor
    function = expected_mass_change01
    execute_on = 'timestep_end'
    outputs = 'console'
  [../]
  [./m01_expect]
    type = FunctionValuePostprocessor
    function = mass01_expect
    execute_on = 'timestep_end'
  [../]
  [./p11]
    type = PointValue
    point = '1 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
[]

[BCs]
  [./flux]
    type = PorousFlowSink
    boundary = 'left'
    variable = pp
    use_mobility = true
    use_relperm = true
    fluid_phase = 0
    flux_function = 6
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-3
  end_time = 0.03
  nl_rel_tol = 1E-12
  nl_abs_tol = 1E-12
[]

[Outputs]
  file_base = s02
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
    interval = 30
  [../]
  [./csv]
    type = CSV
    execute_on = 'timestep_end'
    interval = 3
  [../]
[]
(modules/porous_flow/test/tests/sinks/s02.i)

In this case, the expected result is (6) MOOSE correctly produces this result, as illustrated in Figure 3.

Figure 3: Results of Test 3, illustrating that MOOSE correctly applies a constant sink flux modified by the fluid relative permeability.

Test 4

A similar setup to Test 1 is used here, but with a 3-component, single-phase fluid, with the sink flux only extracting the second component, with a rate proportional to the mass fraction of that component.

The input file:

# apply a sink flux on just one component of a 3-component system and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp frac0 frac1'
    number_fluid_phases = 1
    number_fluid_components = 3
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1.1
  [../]
[]

[Variables]
  [./pp]
  [../]
  [./frac0]
    initial_condition = 0.1
  [../]
  [./frac1]
    initial_condition = 0.6
  [../]
[]

[ICs]
  [./pp]
    type = FunctionIC
    variable = pp
    function = y
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = frac0
  [../]
  [./mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = frac1
  [../]
  [./mass2]
    type = PorousFlowMassTimeDerivative
    fluid_component = 2
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'frac0 frac1'
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '0.2 0 0 0 0.1 0 0 0 0.1'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 0
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
[]

[Functions]
  [./mass1_00]
    type = ParsedFunction
    value = 'frac*vol*por*dens0*exp(pp/bulk)*pow(1+pow(-al*pp,1.0/(1-m)),-m)'
    vars = 'frac  vol  por dens0 pp bulk al m'
    vals = 'f1_00 0.25 0.1 1.1  p00 1.3 1.1 0.5'
  [../]
  [./expected_mass_change1_00]
    type = ParsedFunction
    value = 'frac*fcn*area*dt'
    vars = 'frac fcn area dt'
    vals = 'f1_00 6  0.5  1E-3'
  [../]
  [./mass1_00_expect]
    type = ParsedFunction
    value = 'mass_prev-mass_change'
    vars = 'mass_prev mass_change'
    vals = 'm1_00_prev  del_m1_00'
  [../]
  [./mass1_01]
    type = ParsedFunction
    value = 'frac*vol*por*dens0*exp(pp/bulk)*pow(1+pow(-al*pp,1.0/(1-m)),-m)'
    vars = 'frac  vol  por dens0 pp bulk al m'
    vals = 'f1_01 0.25 0.1 1.1  p01 1.3 1.1 0.5'
  [../]
  [./expected_mass_change1_01]
    type = ParsedFunction
    value = 'frac*fcn*area*dt'
    vars = 'frac fcn area dt'
    vals = 'f1_01 6  0.5  1E-3'
  [../]
  [./mass1_01_expect]
    type = ParsedFunction
    value = 'mass_prev-mass_change'
    vars = 'mass_prev mass_change'
    vals = 'm1_01_prev  del_m1_01'
  [../]
[]

[Postprocessors]
  [./f1_00]
    type = PointValue
    point = '0 0 0'
    variable = frac1
    execute_on = 'initial timestep_end'
  [../]
  [./flux_00]
    type = PointValue
    point = '0 0 0'
    variable = flux_out
    execute_on = 'initial timestep_end'
  [../]
  [./p00]
    type = PointValue
    point = '0 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m1_00]
    type = FunctionValuePostprocessor
    function = mass1_00
    execute_on = 'initial timestep_end'
  [../]
  [./m1_00_prev]
    type = FunctionValuePostprocessor
    function = mass1_00
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./del_m1_00]
    type = FunctionValuePostprocessor
    function = expected_mass_change1_00
    execute_on = 'timestep_end'
    outputs = 'console'
  [../]
  [./m1_00_expect]
    type = FunctionValuePostprocessor
    function = mass1_00_expect
    execute_on = 'timestep_end'
  [../]
  [./f1_01]
    type = PointValue
    point = '0 1 0'
    variable = frac1
    execute_on = 'initial timestep_end'
  [../]
  [./flux_01]
    type = PointValue
    point = '0 1 0'
    variable = flux_out
    execute_on = 'initial timestep_end'
  [../]
  [./p01]
    type = PointValue
    point = '0 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m1_01]
    type = FunctionValuePostprocessor
    function = mass1_01
    execute_on = 'initial timestep_end'
  [../]
  [./m1_01_prev]
    type = FunctionValuePostprocessor
    function = mass1_01
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./del_m1_01]
    type = FunctionValuePostprocessor
    function = expected_mass_change1_01
    execute_on = 'timestep_end'
    outputs = 'console'
  [../]
  [./m1_01_expect]
    type = FunctionValuePostprocessor
    function = mass1_01_expect
    execute_on = 'timestep_end'
  [../]
  [./f1_11]
    type = PointValue
    point = '1 1 0'
    variable = frac1
    execute_on = 'initial timestep_end'
  [../]
  [./flux_11]
    type = PointValue
    point = '1 1 0'
    variable = flux_out
    execute_on = 'initial timestep_end'
  [../]
  [./p11]
    type = PointValue
    point = '1 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
[]

[BCs]
  [./flux]
    type = PorousFlowSink
    boundary = 'left'
    variable = frac1
    use_mobility = false
    use_relperm = false
    mass_fraction_component = 1
    fluid_phase = 0
    flux_function = 6
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-3
  end_time = 0.01
  nl_rel_tol = 1E-12
  nl_abs_tol = 1E-12
[]

[Outputs]
  file_base = s07
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
  [../]
  [./csv]
    type = CSV
    execute_on = 'timestep_end'
  [../]
[]
(modules/porous_flow/test/tests/sinks/s07.i)

This test checks that the flux is correctly implemented (see Figure 4) and that the correct fluid component is being withdrawn from the correct nodes. MOOSE produces the expected result.

Figure 4: Results of Test 4, illustrating that MOOSE correctly applies a sink flux of a particular fluid component proportional to the component's mass fraction.

Test 5

A sink is applied to the left edge () of a 3D mesh. A 3-component, 2-phase fluid is used. Call the two phases "water" and "gas". The porepressures are initialised to and . The mass fractions are initialised to in the water phase, and in the gas phase. The water phase is assumed to have density , and the gas phase . A van-Genuchten capillary relationship is used: (7) with Pa, and . The water relative permeaility is assumed to be Corey type with exponent , and the gas phase has exponent (that is , with ).

The sink flux acts only on the second component. It is multiplied by the relative permeability of the gas phase, and the mass fraction of the second component in the gas phase. This is possibly meaningless physically, but acts as a good test of the PorousFlowSink. In this test the mass fractions remain fixed: there is nothing to induce a change of a component from the water phase to the gas phase since the only Kernels used are mass-conservation Kernels that simply demand mass conservation of each fluid component (summed over each phase).

The input file:

# apply a sink flux on just one component of a 3-component, 2-phase system and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater frac_ph0_c0 pgas'
    number_fluid_phases = 2
    number_fluid_components = 3
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1.1
  [../]
[]

[Variables]
  [./pwater]
  [../]
  [./frac_ph0_c0]
    initial_condition = 0.3
  [../]
  [./pgas]
  [../]
[]

[ICs]
  [./pwater]
    type = FunctionIC
    variable = pwater
    function = y
  [../]
  [./pgas]
    type = FunctionIC
    variable = pgas
    function = y+3
  [../]
[]

[Kernels]
  [./mass_c0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = frac_ph0_c0
  [../]
  [./mass_c1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pwater
  [../]
  [./mass_c2]
    type = PorousFlowMassTimeDerivative
    fluid_component = 2
    variable = pgas
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid0]
      type = SimpleFluidProperties
      bulk_modulus = 2.3
      density0 = 1.5
      thermal_expansion = 0
      viscosity = 2.1
    [../]
    [./simple_fluid1]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'frac_ph0_c0 frac_ph0_c1 frac_ph1_c0 frac_ph1_c1'
  [../]
  [./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 = '0.2 0 0 0 0.1 0 0 0 0.1'
  [../]
  [./relperm0]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 1
    phase = 0
  [../]
  [./relperm1]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 1
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
  [./frac_ph0_c1]
    initial_condition = 0.35
  [../]
  [./frac_ph1_c0]
    initial_condition = 0.1
  [../]
  [./frac_ph1_c1]
    initial_condition = 0.8
  [../]
[]

[Functions]
  [./mass1_00]
    type = ParsedFunction
    value = 'fgas*vol*por*dens0gas*exp(pgas/bulkgas)*(1-pow(1+pow(al*(pgas-pwater),1.0/(1-m)),-m))+fwater*vol*por*dens0water*exp(pwater/bulkwater)*(pow(1+pow(al*(pgas-pwater),1.0/(1-m)),-m))'
    vars = 'vol  por dens0gas pgas    pwater    bulkgas al  m   dens0water bulkwater fgas           fwater'
    vals = '0.25 0.1 1.1      pgas_00 pwater_00 1.3     1.1 0.5 1.5        2.3       frac_ph1_c1_00 frac_ph0_c1_00'
  [../]
  [./expected_mass_change1_00]
    type = ParsedFunction
    value = 'frac*fcn*area*dt*pow(1-pow(1+pow(al*(pgas-pwater),1.0/(1-m)),-m), 2)'
    vars = 'frac           fcn area dt  pgas    pwater    al m'
    vals = 'frac_ph1_c1_00 100 0.5 1E-3 pgas_00 pwater_00 1.1 0.5'
  [../]
  [./mass1_00_expect]
    type = ParsedFunction
    value = 'mass_prev-mass_change'
    vars = 'mass_prev mass_change'
    vals = 'm1_00_prev  del_m1_00'
  [../]
[]

[Postprocessors]
  [./total_mass_comp0]
    type = PorousFlowFluidMass
    fluid_component = 0
  [../]
  [./total_mass_comp1]
    type = PorousFlowFluidMass
    fluid_component = 1
  [../]
  [./total_mass_comp2]
    type = PorousFlowFluidMass
    fluid_component = 2
  [../]
  [./frac_ph1_c1_00]
    type = PointValue
    point = '0 0 0'
    variable = frac_ph1_c1
    execute_on = 'initial timestep_end'
  [../]
  [./frac_ph0_c1_00]
    type = PointValue
    point = '0 0 0'
    variable = frac_ph0_c1
    execute_on = 'initial timestep_end'
  [../]
  [./flux_00]
    type = PointValue
    point = '0 0 0'
    variable = flux_out
    execute_on = 'initial timestep_end'
  [../]
  [./pgas_00]
    type = PointValue
    point = '0 0 0'
    variable = pgas
    execute_on = 'initial timestep_end'
  [../]
  [./pwater_00]
    type = PointValue
    point = '0 0 0'
    variable = pwater
    execute_on = 'initial timestep_end'
  [../]
  [./m1_00]
    type = FunctionValuePostprocessor
    function = mass1_00
    execute_on = 'initial timestep_end'
  [../]
  [./m1_00_prev]
    type = FunctionValuePostprocessor
    function = mass1_00
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./del_m1_00]
    type = FunctionValuePostprocessor
    function = expected_mass_change1_00
    execute_on = 'timestep_end'
    outputs = 'console'
  [../]
  [./m1_00_expect]
    type = FunctionValuePostprocessor
    function = mass1_00_expect
    execute_on = 'timestep_end'
  [../]
[]

[BCs]
  [./flux_ph1_c1]
    type = PorousFlowSink
    boundary = 'left'
    variable = pwater # sink applied to the mass_c1 Kernel
    use_mobility = false
    use_relperm = true
    mass_fraction_component = 1
    fluid_phase = 1
    flux_function = 100
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 100 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-3
  end_time = 0.01
  nl_rel_tol = 1E-12
  nl_abs_tol = 1E-12
[]

[Outputs]
  file_base = s08
  exodus = true
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
  [../]
  [./csv]
    type = CSV
    execute_on = 'timestep_end'
  [../]
[]
(modules/porous_flow/test/tests/sinks/s08.i)

The test checks whether MOOSE is correctly applying the sink flux, and that the fluid-component masses at the nodes respond correctly to the flux. Figure 5 demonstrates that MOOSE produces the expected result.

Figure 5: Results of Test 5, illustrating that in a 2-phase system MOOSE correctly applies a sink flux of a particular fluid component proportional to the component's mass fraction and the relative permeaility of the gas phase.

Piecewise-linear sink

A sink flux of strength (8) (measured in kg.m.s) is applied to the right side () of a 3D mesh. A single-phase, single-component fluid is used, and the porepressure is initialised to (for ). No fluid flow within the element is used, so the masses of fluid at the finite-element nodes behave independently. The fluid is assumed to have density kg.m. The porosity is 0.1.

The input file:

# apply a piecewise-linear sink flux and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  [../]
[]

[Variables]
  [./pp]
  [../]
[]

[ICs]
  [./pp]
    type = FunctionIC
    variable = pp
    function = y+1
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-5 0 0 0 1E-5 0 0 0 1E-5'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 0
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
  [./xval]
  [../]
  [./yval]
  [../]
  [./pt_shift]
    initial_condition = 0.3
  [../]
[]

[ICs]
  [./xval]
    type = FunctionIC
    variable = xval
    function = x
  [../]
  [./yval]
    type = FunctionIC
    variable = yval
    function = y
  [../]
[]

[Functions]
  [./mass10]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p10 1.3'
  [../]
  [./rate10]
    type = ParsedFunction
    value = 'fcn*if(pp>0.8,1,if(pp<0.3,0.5,0.2+pp))'
    vars = 'fcn pp'
    vals = '8   p10'
  [../]
  [./mass10_expect]
    type = ParsedFunction
    value = 'mass_prev-rate*area*dt'
    vars = 'mass_prev rate     area dt'
    vals = 'm10_prev  m10_rate 0.5 1E-3'
  [../]
  [./mass11]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)'
    vars = 'vol por dens0 pp bulk'
    vals = '0.25 0.1 1.1 p11 1.3'
  [../]
  [./rate11]
    type = ParsedFunction
    value = 'fcn*if(pp>0.8,1,if(pp<0.3,0.5,0.2+pp))'
    vars = 'fcn pp'
    vals = '8   p11'
  [../]
  [./mass11_expect]
    type = ParsedFunction
    value = 'mass_prev-rate*area*dt'
    vars = 'mass_prev rate     area dt'
    vals = 'm11_prev  m11_rate 0.5 1E-3'
  [../]
[]

[Postprocessors]
  [./p00]
    type = PointValue
    point = '0 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p10]
    type = PointValue
    point = '1 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m10]
    type = FunctionValuePostprocessor
    function = mass10
    execute_on = 'initial timestep_end'
  [../]
  [./m10_prev]
    type = FunctionValuePostprocessor
    function = mass10
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./m10_rate]
    type = FunctionValuePostprocessor
    function = rate10
    execute_on = 'timestep_end'
  [../]
  [./m10_expect]
    type = FunctionValuePostprocessor
    function = mass10_expect
    execute_on = 'timestep_end'
  [../]
  [./p01]
    type = PointValue
    point = '0 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p11]
    type = PointValue
    point = '1 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m11]
    type = FunctionValuePostprocessor
    function = mass11
    execute_on = 'initial timestep_end'
  [../]
  [./m11_prev]
    type = FunctionValuePostprocessor
    function = mass11
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./m11_rate]
    type = FunctionValuePostprocessor
    function = rate11
    execute_on = 'timestep_end'
  [../]
  [./m11_expect]
    type = FunctionValuePostprocessor
    function = mass11_expect
    execute_on = 'timestep_end'
  [../]
[]

[BCs]
  [./flux]
    type = PorousFlowPiecewiseLinearSink
    boundary = 'right'
    PT_shift = pt_shift
    pt_vals = '0.0 0.5'
    multipliers = '0.5 1'
    variable = pp
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 8
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-3
  end_time = 1E-2
  nl_rel_tol = 1E-12
  nl_abs_tol = 1E-12
[]

[Outputs]
  file_base = s04
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
  [../]
  [./csv]
    type = CSV
    execute_on = 'timestep_end'
  [../]
[]
(modules/porous_flow/test/tests/sinks/s04.i)

Under these conditions, the expected result for the fluid mass at a node on the right side of the mesh is (9) The notation is the same as in previous sections.

The test checks that the mass evolves according to this equation, and that the flux is applied correctly. Figure 6 demonstrates agreement with the expected flux and the MOOSE implementation.

Figure 6: A piecewise-linear sink flux is correctly modelled by MOOSE.

Half-Gaussian sink

A sink flux of strength (10) (measured in kg.m.s) is applied to the right side () of a 3D mesh. This is a half-Gaussian sink with center 0.9Pa, standard deviation 0.5Pa and maximum 6. A single-phase, single-component fluid is used, and the porepressure is initialised to (for ). No fluid flow within the element is used, so the masses of fluid at the finite-element nodes behave independently. The fluid is assumed to have density kg.m. The porosity is 0.1. A van-Genuchten capillary relationship is used: (11) with \,Pa, and .

The input file:

# apply a half-gaussian sink flux and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1.1
  [../]
[]

[Variables]
  [./pp]
  [../]
[]

[ICs]
  [./pp]
    type = FunctionIC
    variable = pp
    function = y+1.4
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-5 0 0 0 1E-5 0 0 0 1E-5'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 0
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
[]

[Functions]
  [./mass10]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)*if(pp>=0,1,pow(1+pow(-al*pp,1.0/(1-m)),-m))'
    vars = 'vol por dens0 pp bulk al m'
    vals = '0.25 0.1 1.1 p10 1.3 1.1 0.5'
  [../]
  [./rate10]
    type = ParsedFunction
    value = 'if(pp>center,fcn,fcn*exp(-0.5*(pp-center)*(pp-center)/sd/sd))'
    vars = 'fcn pp  center sd'
    vals = '6   p10 0.9    0.5'
  [../]
  [./mass10_expect]
    type = ParsedFunction
    value = 'mass_prev-rate*area*dt'
    vars = 'mass_prev rate     area dt'
    vals = 'm10_prev  m10_rate 0.5 2E-3'
  [../]
  [./mass11]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)*if(pp>=0,1,pow(1+pow(-al*pp,1.0/(1-m)),-m))'
    vars = 'vol por dens0 pp bulk al m'
    vals = '0.25 0.1 1.1 p11 1.3 1.1 0.5'
  [../]
  [./rate11]
    type = ParsedFunction
    value = 'if(pp>center,fcn,fcn*exp(-0.5*(pp-center)*(pp-center)/sd/sd))'
    vars = 'fcn pp  center sd'
    vals = '6   p11 0.9    0.5'
  [../]
  [./mass11_expect]
    type = ParsedFunction
    value = 'mass_prev-rate*area*dt'
    vars = 'mass_prev rate     area dt'
    vals = 'm11_prev  m11_rate 0.5 2E-3'
  [../]
[]

[Postprocessors]
  [./flux10]
    type = PointValue
    variable = flux_out
    point = '1 0 0'
  [../]
  [./p00]
    type = PointValue
    point = '0 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p10]
    type = PointValue
    point = '1 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m10]
    type = FunctionValuePostprocessor
    function = mass10
    execute_on = 'initial timestep_end'
  [../]
  [./m10_prev]
    type = FunctionValuePostprocessor
    function = mass10
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./m10_rate]
    type = FunctionValuePostprocessor
    function = rate10
    execute_on = 'timestep_end'
  [../]
  [./m10_expect]
    type = FunctionValuePostprocessor
    function = mass10_expect
    execute_on = 'timestep_end'
  [../]
  [./p01]
    type = PointValue
    point = '0 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p11]
    type = PointValue
    point = '1 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m11]
    type = FunctionValuePostprocessor
    function = mass11
    execute_on = 'initial timestep_end'
  [../]
  [./m11_prev]
    type = FunctionValuePostprocessor
    function = mass11
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./m11_rate]
    type = FunctionValuePostprocessor
    function = rate11
    execute_on = 'timestep_end'
  [../]
  [./m11_expect]
    type = FunctionValuePostprocessor
    function = mass11_expect
    execute_on = 'timestep_end'
  [../]
[]

[BCs]
  [./flux]
    type = PorousFlowHalfGaussianSink
    boundary = 'right'
    max = 6
    sd = 0.5
    center = 0.9
    variable = pp
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

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

[Outputs]
  file_base = s05
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
    interval = 5
  [../]
  [./csv]
    type = CSV
    execute_on = 'timestep_end'
    interval = 3
  [../]
[]
(modules/porous_flow/test/tests/sinks/s05.i)

Under these conditions, the expected result for the fluid mass at a node on the right side of the mesh is (12) The notation is the same as in previous sections.

The test checks that the mass evolves according to this equation, and that the flux is applied correctly. Figure 7 demonstrates agreement with the expected flux and the MOOSE implementation.

Figure 7: A half-Gaussian sink flux with center 0.9Pa and standard deviation 0.5Pa is correctly modelled by MOOSE.

Half-cubic sink

A sink flux of strength (13) (measured in kg.m.s) is applied to the right side () of a 3D mesh. This is a half-cubic sink with center 0.9Pa, cutoff Pa, and maximum 3kg.m.s. A single-phase, single-component fluid is used, and the porepressure is initialised to (for and ). No fluid flow within the element is used, so the masses of fluid at the finite-element nodes behave independently. The fluid is assumed to have density kg.m. The porosity is 0.1.

The input file:

# apply a half-cubic sink flux and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 1
  ny = 1
  nz = 1
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1.1
  [../]
[]

[Variables]
  [./pp]
  [../]
[]

[ICs]
  [./pp]
    type = FunctionIC
    variable = pp
    function = x*(y+1)
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1.3
      density0 = 1.1
      thermal_expansion = 0
      viscosity = 1.1
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-5 0 0 0 1E-5 0 0 0 1E-5'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    at_nodes = true
    n = 2
    phase = 0
  [../]
[]

[AuxVariables]
  [./flux_out]
  [../]
[]

[Functions]
  [./mass10]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)*if(pp>=0,1,pow(1+pow(-al*pp,1.0/(1-m)),-m))'
    vars = 'vol por dens0 pp bulk al m'
    vals = '0.25 0.1 1.1 p10 1.3 1.1 0.5'
  [../]
  [./rate10]
    type = ParsedFunction
    value = 'fcn*if(pp>center,m,if(pp<themin,0,m/c/c/c*(2*(pp-center)+c)*((pp-center)-c)*((pp-center)-c)))'
    vars = 'm fcn pp  center sd  themin c'
    vals = '2 3   p10 0.9    0.5 0.1   -0.8'
  [../]
  [./mass10_expect]
    type = ParsedFunction
    value = 'mass_prev-rate*area*dt'
    vars = 'mass_prev rate     area dt'
    vals = 'm10_prev  m10_rate 0.5 2E-3'
  [../]
  [./mass11]
    type = ParsedFunction
    value = 'vol*por*dens0*exp(pp/bulk)*if(pp>=0,1,pow(1+pow(-al*pp,1.0/(1-m)),-m))'
    vars = 'vol por dens0 pp bulk al m'
    vals = '0.25 0.1 1.1 p11 1.3 1.1 0.5'
  [../]
  [./rate11]
    type = ParsedFunction
    value = 'fcn*if(pp>center,m,if(pp<themin,0,m/c/c/c*(2*(pp-center)+c)*((pp-center)-c)*((pp-center)-c)))'
    vars = 'm fcn pp  center sd  themin c'
    vals = '2 3   p11 0.9    0.5 0.1   -0.8'
  [../]
  [./mass11_expect]
    type = ParsedFunction
    value = 'mass_prev-rate*area*dt'
    vars = 'mass_prev rate     area dt'
    vals = 'm11_prev  m11_rate 0.5 2E-3'
  [../]
[]

[Postprocessors]
  [./flux00]
    type = PointValue
    variable = flux_out
    point = '0 0 0'
  [../]
  [./flux01]
    type = PointValue
    variable = flux_out
    point = '0 1 0'
  [../]
  [./flux10]
    type = PointValue
    variable = flux_out
    point = '1 0 0'
  [../]
  [./flux11]
    type = PointValue
    variable = flux_out
    point = '1 1 0'
  [../]
  [./p00]
    type = PointValue
    point = '0 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p10]
    type = PointValue
    point = '1 0 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m10]
    type = FunctionValuePostprocessor
    function = mass10
    execute_on = 'initial timestep_end'
  [../]
  [./m10_prev]
    type = FunctionValuePostprocessor
    function = mass10
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./m10_rate]
    type = FunctionValuePostprocessor
    function = rate10
    execute_on = 'timestep_end'
  [../]
  [./m10_expect]
    type = FunctionValuePostprocessor
    function = mass10_expect
    execute_on = 'timestep_end'
  [../]
  [./p01]
    type = PointValue
    point = '0 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./p11]
    type = PointValue
    point = '1 1 0'
    variable = pp
    execute_on = 'initial timestep_end'
  [../]
  [./m11]
    type = FunctionValuePostprocessor
    function = mass11
    execute_on = 'initial timestep_end'
  [../]
  [./m11_prev]
    type = FunctionValuePostprocessor
    function = mass11
    execute_on = 'timestep_begin'
    outputs = 'console'
  [../]
  [./m11_rate]
    type = FunctionValuePostprocessor
    function = rate11
    execute_on = 'timestep_end'
  [../]
  [./m11_expect]
    type = FunctionValuePostprocessor
    function = mass11_expect
    execute_on = 'timestep_end'
  [../]
[]

[BCs]
  [./flux]
    type = PorousFlowHalfCubicSink
    boundary = 'left right'
    max = 2
    cutoff = -0.8
    center = 0.9
    variable = pp
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 3
    save_in = flux_out
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

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

[Outputs]
  file_base = s06
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
    interval = 5
  [../]
  [./csv]
    type = CSV
    execute_on = 'timestep_end'
    interval = 3
  [../]
[]
(modules/porous_flow/test/tests/sinks/s06.i)

Under these conditions, the expected result for the fluid mass at a node on the right side of the mesh is (14) The notation is the same as in previous sections.

The test checks that the mass evolves according to this equation, and that the flux is applied correctly. Figure 8 demonstrates agreement with the expected flux and the MOOSE implementation.

Figure 8: A half-cubic sink with center 0.9Pa, cutoff -0.8Pa, and maximum 3kg/m2/s is correctly modelled by MOOSE.

Piecewise linear sink, and advection of a fluid component

The porepressure at a boundary may be maintained at a fixed value by applying a sufficiently strong piecewise-linear sink: (15) (measured in kg.m.s) for large conductance . Note that if is too large then it will dominate the numerics and MOOSE will not converge.

Similarly, for the multi-component case, the flux of of fluid component should be made proportional to the component mass fraction, : (16) This is a "natural" boundary condition, in that fluid exits or enters the porous material at a rate dictated by the mass-fraction within the porous material. This means, for instance, that if fluid is exiting ( in this case) then only components that exist at the boundary system will exit, and MOOSE will not attempt to extract fluid components that have zero mass-fraction.

This example concerns a 1D porous material occupying the space . It contains a single phase fluid with two fluid components. The porous material initially only contains fluid component , and there is a pressure gradient: (17) For , fluid component is introduced on the material's left side (), by applying the fixed boundary conditions: (18) The right-hand side, at , is subjected to the flux of Eq. (16). The fluid-component flows from the left side to the right side via the pressure gradient. To simplify the following analysis, the fluid bulk modulus is taken to be very large.

The input file:

# Apply a piecewise-linear sink flux to the right-hand side and watch fluid flow to it
#
# This test has a single phase with two components.  The test initialises with
# the porous material fully filled with component=1.  The left-hand side is fixed
# at porepressure=1 and mass-fraction of the zeroth component being unity.
# The right-hand side has a very strong piecewise-linear flux that keeps the
# porepressure~0 at that side.  Fluid mass is extracted by this flux in proportion
# to the fluid component mass fraction.
#
# Therefore, the zeroth fluid component will flow from left to right (down the
# pressure gradient).
#
# The important DE is
# porosity * dc/dt = (perm / visc) * grad(P) * grad(c)
# which is true for c = mass-fraction, and very large bulk modulus of the fluid.
# For grad(P) constant in time and space (as in this example) this is just the
# advection equation for c, with velocity = perm / visc / porosity.  The parameters
# are chosen to velocity = 1 m/s.
# In the numerical world, and especially with full upwinding, the advection equation
# suffers from diffusion.  In this example, the diffusion is obvious when plotting
# the mass-fraction along the line, but the average velocity of the front is still
# correct at 1 m/s.
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmin = 0
  xmax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp frac'
    number_fluid_phases = 1
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  [../]
[]

[Variables]
  [./pp]
  [../]
  [./frac]
  [../]
[]

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

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = frac
  [../]
  [./mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pp
  [../]
  [./flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    gravity = '0 0 0'
    variable = frac
  [../]
  [./flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    gravity = '0 0 0'
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1e10 # need large in order for constant-velocity advection
      density0 = 1 # almost irrelevant, except that the ability of the right BC to keep P fixed at zero is related to density_P0
      thermal_expansion = 0
      viscosity = 11
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = frac
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.1 0 0 0 1.1 0 0 0 1.1'
  [../]
  [./relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2 # irrelevant in this fully-saturated situation
    phase = 0
  [../]
[]

[BCs]
  [./lhs_fixed_a]
    type = PresetBC
    boundary = 'left'
    variable = frac
    value = 1
  [../]
  [./lhs_fixed_b]
    type = PresetBC
    boundary = 'left'
    variable = pp
    value = 1
  [../]
  [./flux0]
    type = PorousFlowPiecewiseLinearSink
    boundary = 'right'
    pt_vals = '-100 100'
    multipliers = '-1 1'
    variable = frac # the zeroth comonent
    mass_fraction_component = 0
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1E4
  [../]
  [./flux1]
    type = PorousFlowPiecewiseLinearSink
    boundary = 'right'
    pt_vals = '-100 100'
    multipliers = '-1 1'
    variable = pp # comonent 1
    mass_fraction_component = 1
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1E4
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-2
  end_time = 1
  nl_rel_tol = 1E-12
  nl_abs_tol = 1E-12
[]

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

[Outputs]
  file_base = s09
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
  [../]
  [./csv]
    type = CSV
    sync_times = '0.1 0.5 1'
    sync_only = true
  [../]
  interval = 10
[]
(modules/porous_flow/test/tests/sinks/s09.i)

Because the fluid bulk modulus is very large, is a solution for all time. This means that the governing equation reduces to (19) In this equation is the porosity, is the permeability tensor, and is the fluid viscosity. This is just the advection equation for the mass fraction , with velocity (20) In this test, the parameters are chosen such that velocitym.s.

The sharp front (described by the advection equation with the initial and boundary conditions) is not maintained by MOOSE. This is due to numerical diffusion, which is particularly strong in the upwinding scheme used in this test. For many more details, see the stabilization page. Nevertheless, MOOSE advects the smooth front with the correct velocity, as shown in Figure 9.

The sharp front is not maintained by MOOSE even when no upwinding is used. In the case at hand, which uses a fully-saturated single-phase fluid, the FullySaturated versions of the Kernels may be used in order to compare with the standard fully-upwinded Kernels. These Kernels do not employ any upwinding whatsoever, so less numerical diffusion is expected. This is demonstrated in Figure 9. Two additional points may also be nocied: (1) the lack of upwinding has produced a "bump" in the mass-fraction profile near the concentrated side; (2) the lack of upwinding means the temperature profile moves slightly slower than it should. These two affects reduce as the mesh density is increased, however.

The input file using the fully-saturated Kernels:

# Apply a piecewise-linear sink flux to the right-hand side and watch fluid flow to it
#
# This test has a single phase with two components.  The test initialises with
# the porous material fully filled with component=1.  The left-hand side is fixed
# at porepressure=1 and mass-fraction of the zeroth component being unity.
# The right-hand side has a very strong piecewise-linear flux that keeps the
# porepressure~0 at that side.  Fluid mass is extracted by this flux in proportion
# to the fluid component mass fraction.
#
# Therefore, the zeroth fluid component will flow from left to right (down the
# pressure gradient).
#
# The important DE is
# porosity * dc/dt = (perm / visc) * grad(P) * grad(c)
# which is true for c = mass-fraction, and very large bulk modulus of the fluid.
# For grad(P) constant in time and space (as in this example) this is just the
# advection equation for c, with velocity = perm / visc / porosity.  The parameters
# are chosen to velocity = 1 m/s.
# In the numerical world, and especially with full upwinding, the advection equation
# suffers from diffusion.  In this example, the diffusion is obvious when plotting
# the mass-fraction along the line, but the average velocity of the front is still
# correct at 1 m/s.
# This test uses the FullySaturated version of the flow Kernel.  This does not
# suffer from as much numerical diffusion as the standard PorousFlow Kernel since
# it does not employ any upwinding.
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmin = 0
  xmax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp frac'
    number_fluid_phases = 1
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
  [../]
[]

[Variables]
  [./pp]
  [../]
  [./frac]
  [../]
[]

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

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = frac
  [../]
  [./mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pp
  [../]
  [./flux0]
    type = PorousFlowFullySaturatedDarcyFlow
    fluid_component = 0
    gravity = '0 0 0'
    variable = frac
  [../]
  [./flux1]
    type = PorousFlowFullySaturatedDarcyFlow
    fluid_component = 1
    gravity = '0 0 0'
    variable = pp
  [../]
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 1e10 # need large in order for constant-velocity advection
      density0 = 1 # almost irrelevant, except that the ability of the right BC to keep P fixed at zero is related to density_P0
      thermal_expansion = 0
      viscosity = 11
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
  [../]
  [./ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = frac
  [../]
  [./simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    at_nodes = false
    permeability = '1.1 0 0 0 1.1 0 0 0 1.1'
  [../]
  [./relperm_qp]
    type = PorousFlowRelativePermeabilityCorey
    n = 2 # irrelevant in this fully-saturated situation
    phase = 0
  [../]
[]

[BCs]
  [./lhs_fixed_a]
    type = PresetBC
    boundary = 'left'
    variable = frac
    value = 1
  [../]
  [./lhs_fixed_b]
    type = PresetBC
    boundary = 'left'
    variable = pp
    value = 1
  [../]
  [./flux0]
    type = PorousFlowPiecewiseLinearSink
    boundary = 'right'
    pt_vals = '-100 100'
    multipliers = '-1 1'
    variable = frac # the zeroth comonent
    mass_fraction_component = 0
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1E4
  [../]
  [./flux1]
    type = PorousFlowPiecewiseLinearSink
    boundary = 'right'
    pt_vals = '-100 100'
    multipliers = '-1 1'
    variable = pp # comonent 1
    mass_fraction_component = 1
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 1E4
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E-2
  end_time = 1
  nl_rel_tol = 1E-11
  nl_abs_tol = 1E-11
[]

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

[Outputs]
  file_base = s09_fully_saturated
  [./console]
    type = Console
    execute_on = 'nonlinear linear'
  [../]
  [./csv]
    type = CSV
    sync_times = '0.1 0.5 1'
    sync_only = true
  [../]
  interval = 10
[]
(modules/porous_flow/test/tests/sinks/s09_fully_saturated.i)

Figure 9: Results of the advection of a fluid component test, illustrating that the numerical implementation of porous flow within MOOSE diffuses sharp fronts, but advects them at the correct velocity (which is 1m/s in this case). Notice the centre of the front is at the correct position in each picture. Less diffusion is experienced when upwinding is not used, but notice the slight "bump" in the non-upwinded version at early times.