Restarting from previous simulations

MOOSE provides several options for restarting or recovering from previous simulations. In this example, we demonstrate how this capability can be used in the PorousFlow module by considering a simple two-part problem. The first is to establish the hydrostatic pressure gradient in a reservoir that is due to gravity. This is the equilibrium state of the reservoir, and will be used as the initial condition for porepressure in the second part of this problem, where gas is injected into the reservoir in a single injection well.

Although this problem is simple, it demonstrates how the results of a previous simulation can be used to provide the initial condition for a second simulation.

Model

The properties of the reservoir used in the model are summarized in Table 1. The pressure gradient is representative of a reservoir at a depth of approximately 1,000 m.

Table 1: Reservoir properties

PropertyValue
Height100 m
Length5,000 m
PressureHydrostatic (10 - 11 MPa)
Temperature50 C
Permeability m (~100 md)
Porosity0.1
NaCl mass fraction0.1
Methane injection rate1 kg/s

Relative permeability of the aqueous phase is represented by a Corey model, with exponent and residual saturation . Gas phase relative permeability is also represented using a Corey model, with exponent and residual saturation . Capillary pressure is given by a van Genuchten model, with exponent , and residual saturation .

Gravity equilibrium

The hydrostatic pressure gradient in the reservoir is the result of gravity and the density of the resident brine. This pressure gradient can be calculated using a fully saturated single-phase model with a relatively coarse mesh. As this pressure gradient depends on gravity, it develops in the vertical direction only, and can therefore be modelled using a two-dimensional Cartesian mesh.

The following input file to establish hydrostatic equilibrium in this model is provided:

# Initial run to establish gravity equilibrium. As only brine is present (no gas),
# we can use the single phase equation of state and kernels, reducing the computational
# cost. An estimate of the hydrostatic pressure gradient is used as the initial condition
# using an approximate brine density of 1060 kg/m^3.
# The end time is set to a large value (~100 years) to allow the pressure to reach
# equilibrium. Steady state detection is used to halt the run when a steady state is reached.

[Mesh]
  type = GeneratedMesh
  dim = 2
  ny = 10
  nx = 10
  ymax = 100
  xmax = 5000
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 -9.81 0'
  temperature_unit = Celsius
[]

[Variables]
  [./porepressure]
  [../]
[]

[ICs]
  [./porepressure]
    type = FunctionIC
    function = ppic
    variable = porepressure
  [../]
[]

[Functions]
  [./ppic]
    type = ParsedFunction
    value = '10e6 + 1060*9.81*(100-y)'
  [../]
[]

[BCs]
  [./top]
    type = PresetBC
    variable = porepressure
    value = 10e6
    boundary = top
  [../]
[]

[AuxVariables]
  [./temperature]
    initial_condition = 50
  [../]
  [./xnacl]
    initial_condition = 0.1
  [../]
  [./brine_density]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    variable = porepressure
  [../]
  [./flux0]
    type = PorousFlowFullySaturatedDarcyFlow
    variable = porepressure
  [../]
[]

[AuxKernels]
  [./brine_density]
    type = PorousFlowPropertyAux
    property = density
    variable = brine_density
    execute_on = 'initial timestep_end'
  [../]
[]

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

[Modules]
  [./FluidProperties]
    [./brine]
      type = BrineFluidProperties
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
    temperature = temperature
  [../]
  [./ps]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
  [../]
  [./brine]
    type = PorousFlowBrine
    compute_enthalpy = false
    compute_internal_energy = false
    xnacl = xnacl
    phase = 0
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0 0 1e-13 0  0 0 1e-13'
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 3e9
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-06
  steady_state_detection = true
  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e1
  [../]
[]

[Outputs]
  execute_on = 'initial timestep_end'
  exodus = true
  perf_graph = true
[]
(modules/porous_flow/examples/restart/gravityeq.i)

In this example, an initial guess at the hydrostatic equilibrium is provided using a function with an estimate of the brine density at the given pressure and temperature (approximately 1060 kg/m)

[Functions]
  [./ppic]
    type = ParsedFunction
    value = '10e6 + 1060*9.81*(100-y)'
  [../]
[]
(modules/porous_flow/examples/restart/gravityeq.i)

The simulation is set to run for approximately 100 years, with steady-state detection enabled to end the run when a steady-state solution has been found

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 3e9
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-06
  steady_state_detection = true
  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e1
  [../]
[]
(modules/porous_flow/examples/restart/gravityeq.i)

This input file quickly establishes gravity equilibrium in twelve time steps, with a hydrostatic pressure gradient shown in Figure 1.

Figure 1: Hydrostatic pressure gradient at gravity equilibrium

Basic restart

The results of the gravity equilibrium simulation can be used as initial conditions for a more complicated simulation. In this case, we now consider injection of methane through an injection well in a two-dimensional radial model constructed by rotating the mesh used in the gravity equilibrium run around the Y-axis.

In contrast to the simple gravity equilibrium model, we now use a two-phase model, with the initial condition for the liquid porepressure taken from the final results of the gravity equilibrium simulation.

As the mesh is identical in both cases (despite the rotation about the Y-axis), we can use the simple variable initialisation available in MOOSE.

The input file for this case is

# Using the results from the equilibrium run to provide the initial condition for
# porepressure, we now inject a gas phase into the brine-saturated reservoir. In this
# example, where the mesh used is identical to the mesh used in gravityeq.i, we can use
# the basic restart capability by simply setting the initial condition for porepressure
# using the results from gravityeq.i.
#
# Even though the gravity equilibrium is established using a 2D mesh, in this example,
# we shift the mesh 0.1 m to the right and rotate it about the Y axis to make a 2D radial
# model.
#
# Methane injection takes place over the surface of the hole created by rotating the mesh,
# and hence the injection area is 2 pi r h. We can calculate this using an AreaPostprocessor,
# and then use this in a ParsedFunction to calculate the injection rate so that 10 kg/s of
# methane is injected.
#
# Results can be improved by uniformly refining the initial mesh.
#
# Note: as this example uses the results from a previous simulation, gravityeq.i MUST be
# run before running this input file.

[Mesh]
  type = FileMesh
  file = gravityeq_out.e
  uniform_refine = 1
[]

[MeshModifiers]
  [./translate]
    type = Transform
    transform = TRANSLATE
    vector_value = '0.1 0 0'
  [../]
[]

[Problem]
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 -9.81 0'
  temperature_unit = Celsius
[]

[Variables]
  [./pp_liq]
    initial_from_file_var = porepressure
  [../]
  [./sat_gas]
    initial_condition = 0
  [../]
[]

[AuxVariables]
  [./temperature]
    initial_condition = 50
  [../]
  [./xnacl]
    initial_condition = 0.1
  [../]
  [./brine_density]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./methane_density]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./massfrac_ph0_sp0]
    initial_condition = 1
  [../]
  [./massfrac_ph1_sp0]
    initial_condition = 0
  [../]
  [./pp_gas]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./sat_liq]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    variable = pp_liq
  [../]
  [./flux0]
    type = PorousFlowAdvectiveFlux
    variable = pp_liq
  [../]
  [./mass1]
    type = PorousFlowMassTimeDerivative
    variable = sat_gas
    fluid_component = 1
  [../]
  [./flux1]
    type = PorousFlowAdvectiveFlux
    variable = sat_gas
    fluid_component = 1
  [../]
[]

[AuxKernels]
  [./brine_density]
    type = PorousFlowPropertyAux
    property = density
    variable = brine_density
    execute_on = 'initial timestep_end'
  [../]
  [./methane_density]
    type = PorousFlowPropertyAux
    property = density
    variable = methane_density
    phase = 1
    execute_on = 'initial timestep_end'
  [../]
  [./pp_gas]
    type = PorousFlowPropertyAux
    property = pressure
    phase = 1
    variable = pp_gas
    execute_on = 'initial timestep_end'
  [../]
  [./sat_liq]
    type = PorousFlowPropertyAux
    property = saturation
    variable = sat_liq
    execute_on = 'initial timestep_end'
  [../]
[]

[BCs]
  [./gas_injection]
    type = PorousFlowSink
    boundary = left
    variable = sat_gas
    flux_function = injection_rate
    fluid_phase = 1
  [../]
  [./brine_out]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pp_liq
    multipliers = '0 1e9'
    pt_vals = '0 1e9'
    fluid_phase = 0
    flux_function = 1e-6
    use_mobility = true
  [../]
[]

[Functions]
  [./injection_rate]
    type = ParsedFunction
    vals = injection_area
    vars = area
    value = '-10/area'
  [../]
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp_liq sat_gas'
    number_fluid_phases = 2
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1e-5
    m = 0.5
    sat_lr = 0.2
  [../]
[]

[Modules]
  [./FluidProperties]
    [./brine]
      type = BrineFluidProperties
    [../]
    [./methane]
      type = MethaneFluidProperties
    [../]
    [./methane_tab]
      type = TabulatedFluidProperties
      fp = methane
      save_file = false
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
    temperature = temperature
  [../]
  [./ps]
    type = PorousFlow2PhasePS
    phase0_porepressure = pp_liq
    phase1_saturation = sat_gas
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  [../]
  [./brine]
    type = PorousFlowBrine
    compute_enthalpy = false
    compute_internal_energy = false
    xnacl = xnacl
    phase = 0
  [../]
  [./methane]
    type = PorousFlowSingleComponentFluid
    compute_enthalpy = false
    compute_internal_energy = false
    fp = methane_tab
    phase = 1
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0 0 1e-13 0  0 0 1e-13'
  [../]
  [./relperm_liq]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.2
    sum_s_res = 0.3
  [../]
  [./relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
    s_res = 0.1
    sum_s_res = 0.3
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
    petsc_options_value = ' asm      lu           NONZERO'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1e8
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-06
  nl_max_its = 20
  dtmax = 1e6
  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e1
  [../]
[]

[Postprocessors]
  [./mass_ph0]
    type = PorousFlowFluidMass
    fluid_component = 0
    execute_on = 'initial timestep_end'
  [../]
  [./mass_ph1]
    type = PorousFlowFluidMass
    fluid_component = 1
    execute_on = 'initial timestep_end'
  [../]
  [./injection_area]
    type = AreaPostprocessor
    boundary = left
    execute_on = initial
  [../]
[]

[Outputs]
  execute_on = 'initial timestep_end'
  exodus = true
  perf_graph = true
[]
(modules/porous_flow/examples/restart/gas_injection.i)

The input mesh for this problem is the output mesh from the gravity equilibrium simulation, which we uniformly refine once to increase the spatial resolution, translate horizontally and rotate about the Y-axis to form a two-dimensional mesh with an injection well at the centre

[Mesh]
  type = FileMesh
  file = gravityeq_out.e
  uniform_refine = 1
[]

[MeshModifiers]
  [./translate]
    type = Transform
    transform = TRANSLATE
    vector_value = '0.1 0 0'
  [../]
[]

[Problem]
  coord_type = RZ
  rz_coord_axis = Y
[]
(modules/porous_flow/examples/restart/gas_injection.i)

The initial condition for the liquid porepressure is read directly from the input mesh

[Variables]
  [./pp_liq]
    initial_from_file_var = porepressure
  [../]
  [./sat_gas]
    initial_condition = 0
  [../]
[]
(modules/porous_flow/examples/restart/gas_injection.i)

Running this model for a total of 10 s, we obtain the gas saturation profile shown in Figure 2.

Figure 2: Gas saturation after 10 seconds

Due to the large density contrast between the injection methane (with a density of 175 kg m compared to approximately 1060 kg m), the methane migrates vertically and pools beneath the top of the reservoir, extending to a radial distance of approximately 900 m. Methane immobilised at residual saturation can be observed in the lower part of the reservoir (where the gas saturation is approximately 0.1).

Flexible restart using a SolutionUserObject

In the above example, the mesh used in the gas injection problem is initially identical to the mesh used in the gravity equilibrium run. MOOSE provides a more flexible restart capability where the mesh used in the subsequent simulation does not have to be identical to the mesh used in the initial simulation. Instead, the variable values from the initial simulation can be projected onto a new mesh in a subsequent simulation using a SolutionUserObject and SolutionFunction.

It is clear from the results shown above that the methane is expected near the upper left of the model (near the top of the injection well). We could refine the original mesh in that region using mesh adaptivity, but we instead construct a new mesh that is refined in that region and use a SolutionUserObject and SolutionFunction to project the initial value for liquid porepressure from the coarse gravity equilibrium mesh to the new mesh.

The input file for this example is

# Using the results from the equilibrium run to provide the initial condition for
# porepressure, we now inject a gas phase into the brine-saturated reservoir. In this
# example, the mesh is not identical to the mesh used in gravityeq.i. Rather, it is
# generated so that it is more refined near the injection boundary and at the top of
# the model, as that is where the gas plume will be present.
#
# To use the hydrostatic pressure calculated using the gravity equilibrium run as the initial
# condition for the pressure, a SolutionUserObject is used, along with a SolutionFunction to
# interpolate the pressure from the gravity equilibrium run to the initial condition for liqiud
# porepressure in this example.
#
# Even though the gravity equilibrium is established using a 2D mesh, in this example,
# we use a mesh shifted 0.1 m to the right and rotate it about the Y axis to make a 2D radial
# model.
#
# Methane injection takes place over the surface of the hole created by rotating the mesh,
# and hence the injection area is 2 pi r h. We can calculate this using an AreaPostprocessor,
# and then use this in a ParsedFunction to calculate the injection rate so that 10 kg/s of
# methane is injected.
#
# Note: as this example uses the results from a previous simulation, gravityeq.i MUST be
# run before running this input file.

[Mesh]
  type = GeneratedMesh
  dim = 2
  ny = 25
  nx = 50
  ymax = 100
  xmin = 0.1
  xmax = 5000
  bias_x = 1.05
  bias_y = 0.95
[]

[Problem]
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 -9.81 0'
  temperature_unit = Celsius
[]

[Variables]
  [./pp_liq]
  [../]
  [./sat_gas]
    initial_condition = 0
  [../]
[]

[ICs]
  [./ppliq_ic]
    type = FunctionIC
    variable = pp_liq
    function = ppliq_ic
  [../]
[]

[AuxVariables]
  [./temperature]
    initial_condition = 50
  [../]
  [./xnacl]
    initial_condition = 0.1
  [../]
  [./brine_density]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./methane_density]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./massfrac_ph0_sp0]
    initial_condition = 1
  [../]
  [./massfrac_ph1_sp0]
    initial_condition = 0
  [../]
  [./pp_gas]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./sat_liq]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    variable = pp_liq
  [../]
  [./flux0]
    type = PorousFlowAdvectiveFlux
    variable = pp_liq
  [../]
  [./mass1]
    type = PorousFlowMassTimeDerivative
    variable = sat_gas
    fluid_component = 1
  [../]
  [./flux1]
    type = PorousFlowAdvectiveFlux
    variable = sat_gas
    fluid_component = 1
  [../]
[]

[AuxKernels]
  [./brine_density]
    type = PorousFlowPropertyAux
    property = density
    variable = brine_density
    execute_on = 'initial timestep_end'
  [../]
  [./methane_density]
    type = PorousFlowPropertyAux
    property = density
    variable = methane_density
    phase = 1
    execute_on = 'initial timestep_end'
  [../]
  [./pp_gas]
    type = PorousFlowPropertyAux
    property = pressure
    phase = 1
    variable = pp_gas
    execute_on = 'initial timestep_end'
  [../]
  [./sat_liq]
    type = PorousFlowPropertyAux
    property = saturation
    variable = sat_liq
    execute_on = 'initial timestep_end'
  [../]
[]

[BCs]
  [./gas_injection]
    type = PorousFlowSink
    boundary = left
    variable = sat_gas
    flux_function = injection_rate
    fluid_phase = 1
  [../]
  [./brine_out]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pp_liq
    multipliers = '0 1e9'
    pt_vals = '0 1e9'
    fluid_phase = 0
    flux_function = 1e-6
    use_mobility = true
    use_relperm = true
    mass_fraction_component = 0
  [../]
[]

[Functions]
  [./injection_rate]
    type = ParsedFunction
    vals = injection_area
    vars = area
    value = '-1/area'
  [../]
  [./ppliq_ic]
    type = SolutionFunction
    solution = soln
  [../]
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp_liq sat_gas'
    number_fluid_phases = 2
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1e-5
    m = 0.5
    sat_lr = 0.2
    pc_max = 1e7
  [../]
  [./soln]
    type = SolutionUserObject
    mesh = gravityeq_out.e
    system_variables = porepressure
  [../]
[]

[Modules]
  [./FluidProperties]
    [./brine]
      type = BrineFluidProperties
    [../]
    [./methane]
      type = MethaneFluidProperties
    [../]
    [./methane_tab]
      type = TabulatedFluidProperties
      fp = methane
      save_file = false
    [../]
  [../]
[]

[Materials]
  [./temperature]
    type = PorousFlowTemperature
    temperature = temperature
  [../]
  [./ps]
    type = PorousFlow2PhasePS
    phase0_porepressure = pp_liq
    phase1_saturation = sat_gas
    capillary_pressure = pc
  [../]
  [./massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
  [../]
  [./brine]
    type = PorousFlowBrine
    compute_enthalpy = false
    compute_internal_energy = false
    xnacl = xnacl
    phase = 0
  [../]
  [./methane]
    type = PorousFlowSingleComponentFluid
    compute_enthalpy = false
    compute_internal_energy = false
    fp = methane_tab
    phase = 1
  [../]
  [./porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  [../]
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0 0 5e-14 0  0 0 1e-13'
  [../]
  [./relperm_liq]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.2
    sum_s_res = 0.3
  [../]
  [./relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
    s_res = 0.1
    sum_s_res = 0.3
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
    petsc_options_value = ' asm      lu           NONZERO'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1e8
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-06
  nl_max_its = 20
  dtmax = 1e6
  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e1
    growth_factor = 1.5
  [../]
[]

[Postprocessors]
  [./mass_ph0]
    type = PorousFlowFluidMass
    fluid_component = 0
    execute_on = 'initial timestep_end'
  [../]
  [./mass_ph1]
    type = PorousFlowFluidMass
    fluid_component = 1
    execute_on = 'initial timestep_end'
  [../]
  [./injection_area]
    type = AreaPostprocessor
    boundary = left
    execute_on = initial
  [../]
[]

[Outputs]
  execute_on = 'initial timestep_end'
  exodus = true
  perf_graph = true
[]
(modules/porous_flow/examples/restart/gas_injection_new_mesh.i)

The mesh in this example is biased so that there are more elements in the upper left region, and fewer away from the well

[Mesh]
  type = GeneratedMesh
  dim = 2
  ny = 25
  nx = 50
  ymax = 100
  xmin = 0.1
  xmax = 5000
  bias_x = 1.05
  bias_y = 0.95
[]
(modules/porous_flow/examples/restart/gas_injection_new_mesh.i)

The results from the gravity equilibrium simulation are read using a SolutionUserObject

[./soln]
  type = SolutionUserObject
  mesh = gravityeq_out.e
  system_variables = porepressure
[../]
(modules/porous_flow/examples/restart/gas_injection_new_mesh.i)

A SolutionFunction is used to spatially interpolate these values

[./ppliq_ic]
  type = SolutionFunction
  solution = soln
[../]
(modules/porous_flow/examples/restart/gas_injection_new_mesh.i)

and finally a FunctionIC is used to set the initial condition of the liquid porepressure

[ICs]
  [./ppliq_ic]
    type = FunctionIC
    variable = pp_liq
    function = ppliq_ic
  [../]
[]
(modules/porous_flow/examples/restart/gas_injection_new_mesh.i)

This more refined model can then be run, resulting in the gas saturation profile shown in Figure 3.

Figure 3: Gas saturation after 10 seconds

Due to the additional refinement near the well and the top of the reservoir, the gas saturation profile is now smoother near the well, with much less lateral migration in the lower part of the reservoir and hence much less methane immobilised at residual saturation. As a result, more methane is able to migrate to the top of the reservoir, and the lateral plume extends nearly 200 m further in this model than the coarser model shown above.

Generalisation to more complicated models

The above example of the restart capability demonstrates how PorousFlow can use the results of previous simulations as initial conditions of new simulations, even when changing from single-phase models to a more complex two-phase model. Of course, these approaches can be extended to models where heat and geomechanics are involved by setting initial conditions for additional variables using an identical approach to that demonstrated in this simple example.