Start | Previous | Next

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]
    type = PorousFlowAdvectiveFlux
    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!

  • You can ask on the moose-users google group.

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]
    type = ParsedAddSideset
    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]
    type = PorousFlowAdvectiveFlux
    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]
    type = ParsedAddSideset
    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]
    type = PorousFlowAdvectiveFlux
    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]
    type = ParsedAddSideset
    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]
    type = PorousFlowAdvectiveFlux
    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]
    type = ParsedAddSideset
    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]
    type = PorousFlowAdvectiveFlux
    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!

Start | Previous | Next