# Porous Flow Tutorial Page 10. Unleashing the full power of PorousFlow: using Kernels and Materials

Now we're ready to build an input file from scratch using Kernels and Materials instead of the Actions we've been using so far. Inevitably you'll have to do this yourself when running PorousFlow simulations: you'll always want some PorousFlow features (boundary conditions, line sources, postprocessors, AuxKernels, etc) that aren't built by the Actions.

We're going to build the simulation from Page 08 which involved unsaturated flow. After building the input file, you can compare it with the Action version. The mesh and Variables remain unchanged. For clarity later on, the porepressure variable has been renamed 'pp'

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

(modules/porous_flow/examples/tutorial/10.i)

Let's start by building the PorousFlowDictator. We have one phase, one component and no chemistry. The single PorousFlow variable is called pp:

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = pp
number_fluid_phases = 1
number_fluid_components = 1
[../]

(modules/porous_flow/examples/tutorial/10.i)

It is always useful to put the name of the PorousFlowDictator in the GlobalParams block, since all PorousFlow objects need it:

[GlobalParams]
PorousFlowDictator = dictator
[]

(modules/porous_flow/examples/tutorial/10.i)

The DE is unsaturated single-phase flow. Refer to the governing equations document. Unsaturated single-phase flow is described by the first equation without the coupling terms to solid mechanics, radioactive decay, chemistry and sources:

(1)

Here since there is just one fluid component. The fluid mass is and the fluid velocity is . Hence we have two Kernels (refer to the bottom of governing equations for their types)

[Kernels]
[./time_derivative]
type = PorousFlowMassTimeDerivative
variable = pp
[../]
[./flux]
variable = pp
gravity = '0 0 0'
[../]
[]

(modules/porous_flow/examples/tutorial/10.i)

It is often common to define some AuxVariables for visualisation purposes (or to feed as coupled variables to other MOOSE objects). A useful variable here is the fluid saturation, which is computed using a PorousFlowPropertyAux so the variable must be a constant monomial:

[AuxVariables]
[./sat]
family = MONOMIAL
order = CONSTANT
[../]
[]

[AuxKernels]
[./saturation]
type = PorousFlowPropertyAux
variable = sat
property = saturation
[../]
[]

(modules/porous_flow/examples/tutorial/10.i)

The BCs remain the same: they used a PorousFlowSink to withdraw fluid from the injection area. The fluid properties also remain the same. For reference, these two blocks are:

[BCs]
[./production]
type = PorousFlowSink
variable = pp
fluid_phase = 0
flux_function = 1E-2
use_relperm = true
boundary = injection_area
[../]
[]

[Modules]
[./FluidProperties]
[./the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[../]
[../]
[]

(modules/porous_flow/examples/tutorial/10.i)

The final block to create is the Materials. This is always the most complicated part of creating an input file, and really only comes by experience.

• You can run MOOSE multiple times, each time getting "undefined" errors like the one at the top of Page 09 and slowly add the required Materials until you get no errors. The information near the bottom of Page 09 is useful here.

• You can inspect other similar input files in the PorousFlow test suite. There are over 300 tests so you're very likely to find what you need. Don't expect the files to always correspond to physically-sensible models (sometimes they're deliberately unsensible to test obscure aspects of PorousFlow) but at least they are guaranteed to run!

In this case, the DEs involve porosity, fluid saturation, fluid density, permeability and viscosity. The fluid mass is lumped to the nodes, so we'll only need porosity at the nodes:

# Unsaturated Darcy-Richards flow without using an Action
[Mesh]
type = AnnularMesh
dim = 2
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
tmin = 0
tmax = 1.57079632679
[]

[MeshModifiers]
[./make3D]
type = MeshExtruder
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
[../]
[./shift_down]
type = Transform
transform = TRANSLATE
vector_value = '0 0 -6'
depends_on = make3D
[../]
[./aquifer]
type = SubdomainBoundingBox
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
depends_on = shift_down
[../]
[./injection_area]
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomain_ids = 1
new_sideset_name = 'injection_area'
depends_on = 'aquifer'
[../]
[./rename]
type = RenameBlock
old_block_id = '0 1'
new_block_name = 'caps aquifer'
depends_on = 'injection_area'
[../]
[]

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = pp
number_fluid_phases = 1
number_fluid_components = 1
[../]
[./pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[../]
[]

[GlobalParams]
PorousFlowDictator = dictator
[]

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

[Kernels]
[./time_derivative]
type = PorousFlowMassTimeDerivative
variable = pp
[../]
[./flux]
variable = pp
gravity = '0 0 0'
[../]
[]

[AuxVariables]
[./sat]
family = MONOMIAL
order = CONSTANT
[../]
[]

[AuxKernels]
[./saturation]
type = PorousFlowPropertyAux
variable = sat
property = saturation
[../]
[]

[BCs]
[./production]
type = PorousFlowSink
variable = pp
fluid_phase = 0
flux_function = 1E-2
use_relperm = true
boundary = injection_area
[../]
[]

[Modules]
[./FluidProperties]
[./the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[../]
[../]
[]

[Materials]
[./porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[../]

(modules/porous_flow/examples/tutorial/10.i)

The permeability is also needed:

  [./permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
[../]
[./saturation_calculator]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[../]
[./temperature]
type = PorousFlowTemperature
temperature = 293
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 3
s_res = 0.1
sum_s_res = 0.1
phase = 0
[../]
[]

[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'
[../]
[]

[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-7
[]

[Outputs]
exodus = true
[]

(modules/porous_flow/examples/tutorial/10.i)

These Materials were also in the model of Page 08.

Now we need to build the saturation. This is not computed by the PorousFlowPropertyAux above. That AuxKernel just retrieves the material property and puts it into an AuxVariable. As mentioned in Page 09 PorousFlow simulations can use a variety of primary Variables. A set of fundamental Materials computes porepressures, saturations, temperatures and mass fractions from these primary Variables. In this case, our primary Variable is just porepressure and we only need to compute a single saturation. However, various other PorousFlow objects need temperature as an input, so we must compute it too, even though we just set it to a constant 293. Also, mass fraction needs to be computed, although it is trivially 1.0 for this problem.

# Unsaturated Darcy-Richards flow without using an Action
[Mesh]
type = AnnularMesh
dim = 2
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
tmin = 0
tmax = 1.57079632679
[]

[MeshModifiers]
[./make3D]
type = MeshExtruder
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
[../]
[./shift_down]
type = Transform
transform = TRANSLATE
vector_value = '0 0 -6'
depends_on = make3D
[../]
[./aquifer]
type = SubdomainBoundingBox
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
depends_on = shift_down
[../]
[./injection_area]
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomain_ids = 1
new_sideset_name = 'injection_area'
depends_on = 'aquifer'
[../]
[./rename]
type = RenameBlock
old_block_id = '0 1'
new_block_name = 'caps aquifer'
depends_on = 'injection_area'
[../]
[]

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = pp
number_fluid_phases = 1
number_fluid_components = 1
[../]
[./pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[../]
[]

[GlobalParams]
PorousFlowDictator = dictator
[]

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

[Kernels]
[./time_derivative]
type = PorousFlowMassTimeDerivative
variable = pp
[../]
[./flux]
variable = pp
gravity = '0 0 0'
[../]
[]

[AuxVariables]
[./sat]
family = MONOMIAL
order = CONSTANT
[../]
[]

[AuxKernels]
[./saturation]
type = PorousFlowPropertyAux
variable = sat
property = saturation
[../]
[]

[BCs]
[./production]
type = PorousFlowSink
variable = pp
fluid_phase = 0
flux_function = 1E-2
use_relperm = true
boundary = injection_area
[../]
[]

[Modules]
[./FluidProperties]
[./the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[../]
[../]
[]

[Materials]
[./porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[../]
[./permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
[../]
[./saturation_calculator]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[../]
[./temperature]
type = PorousFlowTemperature
temperature = 293
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 3
s_res = 0.1
sum_s_res = 0.1
phase = 0
[../]
[]

[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'
[../]
[]

[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-7
[]

[Outputs]
exodus = true
[]

(modules/porous_flow/examples/tutorial/10.i)

These are our so-called "fundamental Materials". In almost every PorousFlow simulation you will see similar Materials. The saturation calculator uses the capillary-pressure function. It is a UserObject:

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = pp
number_fluid_phases = 1
number_fluid_components = 1
[../]
[./pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[../]
[]

(modules/porous_flow/examples/tutorial/10.i)

The density and viscosity need to be computed. This is achieved by a PorousFlowSingleComponentFluid

# Unsaturated Darcy-Richards flow without using an Action
[Mesh]
type = AnnularMesh
dim = 2
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
tmin = 0
tmax = 1.57079632679
[]

[MeshModifiers]
[./make3D]
type = MeshExtruder
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
[../]
[./shift_down]
type = Transform
transform = TRANSLATE
vector_value = '0 0 -6'
depends_on = make3D
[../]
[./aquifer]
type = SubdomainBoundingBox
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
depends_on = shift_down
[../]
[./injection_area]
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomain_ids = 1
new_sideset_name = 'injection_area'
depends_on = 'aquifer'
[../]
[./rename]
type = RenameBlock
old_block_id = '0 1'
new_block_name = 'caps aquifer'
depends_on = 'injection_area'
[../]
[]

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = pp
number_fluid_phases = 1
number_fluid_components = 1
[../]
[./pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[../]
[]

[GlobalParams]
PorousFlowDictator = dictator
[]

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

[Kernels]
[./time_derivative]
type = PorousFlowMassTimeDerivative
variable = pp
[../]
[./flux]
variable = pp
gravity = '0 0 0'
[../]
[]

[AuxVariables]
[./sat]
family = MONOMIAL
order = CONSTANT
[../]
[]

[AuxKernels]
[./saturation]
type = PorousFlowPropertyAux
variable = sat
property = saturation
[../]
[]

[BCs]
[./production]
type = PorousFlowSink
variable = pp
fluid_phase = 0
flux_function = 1E-2
use_relperm = true
boundary = injection_area
[../]
[]

[Modules]
[./FluidProperties]
[./the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[../]
[../]
[]

[Materials]
[./porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[../]
[./permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
[../]
[./saturation_calculator]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[../]
[./temperature]
type = PorousFlowTemperature
temperature = 293
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 3
s_res = 0.1
sum_s_res = 0.1
phase = 0
[../]
[]

[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'
[../]
[]

[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-7
[]

[Outputs]
exodus = true
[]

(modules/porous_flow/examples/tutorial/10.i)

which must be joined to a std::vector (see Page 09) with some PorousFlowJoiners:

# Unsaturated Darcy-Richards flow without using an Action
[Mesh]
type = AnnularMesh
dim = 2
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
tmin = 0
tmax = 1.57079632679
[]

[MeshModifiers]
[./make3D]
type = MeshExtruder
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
[../]
[./shift_down]
type = Transform
transform = TRANSLATE
vector_value = '0 0 -6'
depends_on = make3D
[../]
[./aquifer]
type = SubdomainBoundingBox
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
depends_on = shift_down
[../]
[./injection_area]
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomain_ids = 1
new_sideset_name = 'injection_area'
depends_on = 'aquifer'
[../]
[./rename]
type = RenameBlock
old_block_id = '0 1'
new_block_name = 'caps aquifer'
depends_on = 'injection_area'
[../]
[]

[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = pp
number_fluid_phases = 1
number_fluid_components = 1
[../]
[./pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[../]
[]

[GlobalParams]
PorousFlowDictator = dictator
[]

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

[Kernels]
[./time_derivative]
type = PorousFlowMassTimeDerivative
variable = pp
[../]
[./flux]
variable = pp
gravity = '0 0 0'
[../]
[]

[AuxVariables]
[./sat]
family = MONOMIAL
order = CONSTANT
[../]
[]

[AuxKernels]
[./saturation]
type = PorousFlowPropertyAux
variable = sat
property = saturation
[../]
[]

[BCs]
[./production]
type = PorousFlowSink
variable = pp
fluid_phase = 0
flux_function = 1E-2
use_relperm = true
boundary = injection_area
[../]
[]

[Modules]
[./FluidProperties]
[./the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[../]
[../]
[]

[Materials]
[./porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[../]
[./permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
[../]
[./saturation_calculator]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[../]
[./temperature]
type = PorousFlowTemperature
temperature = 293
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 3
s_res = 0.1
sum_s_res = 0.1
phase = 0
[../]
[]

[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'
[../]
[]

[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-7
[]

[Outputs]
exodus = true
[]

(modules/porous_flow/examples/tutorial/10.i)

Finally, the relative permeability must be defined for each phase (just one phase in this example) and Joined into a std::vector

# Unsaturated Darcy-Richards flow without using an Action
[Mesh]
type = AnnularMesh
dim = 2
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
tmin = 0
tmax = 1.57079632679

(modules/porous_flow/examples/tutorial/10.i)

Our input file has been built! You may check that it gives exactly the same answers as the one on Page 08. It is 232 lines long while the one on Page 08 is only 139 lines long. Now that you've reached this point, you can start to build much more powerful models that don't use the PorousFlow Actions: models that involve multi-phase, multi-component flows, complicated thermal couplings, elasto-plasticity, chemical reactions, line sources and sinks, and sophisticated boundary terms!