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
[./TimeStepper]
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
[./TimeStepper]
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]
variable = pp_liq
[../]
[./mass1]
type = PorousFlowMassTimeDerivative
variable = sat_gas
fluid_component = 1
[../]
[./flux1]
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]
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]
variable = pp_liq
[../]
[./mass1]
type = PorousFlowMassTimeDerivative
variable = sat_gas
fluid_component = 1
[../]
[./flux1]
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]
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.