- allow_fp_and_tabulationFalseWhether to allow the two sources of data concurrently
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether to allow the two sources of data concurrently
- create_pT_interpolationsTrueWhether to load (from file) or create (from a fluid property object) properties interpolations from pressure and temperature
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether to load (from file) or create (from a fluid property object) properties interpolations from pressure and temperature
- create_ve_interpolationsFalseWhether to load (from file) or create (from a fluid property object) properties interpolations from pressure and temperature
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether to load (from file) or create (from a fluid property object) properties interpolations from pressure and temperature
- execute_onTIMESTEP_ENDThe list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Unit:(no unit assumed)
Controllable:No
Description:The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
- fluid_property_output_fileName of the CSV file which can be output with the tabulation. This file can then be read as a 'fluid_property_file'
C++ Type:FileName
Unit:(no unit assumed)
Controllable:No
Description:Name of the CSV file which can be output with the tabulation. This file can then be read as a 'fluid_property_file'
- fluid_property_ve_fileName of the csv file containing the tabulated (v,e) fluid property data.
C++ Type:FileName
Unit:(no unit assumed)
Controllable:No
Description:Name of the csv file containing the tabulated (v,e) fluid property data.
- fluid_property_ve_output_fileName of the CSV file which can be output with the (v,e) tabulation. This file can then be read as a 'fluid_property_ve_file'
C++ Type:FileName
Unit:(no unit assumed)
Controllable:No
Description:Name of the CSV file which can be output with the (v,e) tabulation. This file can then be read as a 'fluid_property_ve_file'
- fpThe name of the FluidProperties UserObject
C++ Type:UserObjectName
Unit:(no unit assumed)
Controllable:No
Description:The name of the FluidProperties UserObject
- interpolated_propertiesdensity enthalpy internal_energy viscosityProperties to interpolate if no data file is provided
Default:density enthalpy internal_energy viscosity
C++ Type:MultiMooseEnum
Unit:(no unit assumed)
Controllable:No
Description:Properties to interpolate if no data file is provided
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
- use_log_grid_eFalseOption to use a base-10 logarithmically-spaced grid for specific internal energy instead of a linearly-spaced grid.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Option to use a base-10 logarithmically-spaced grid for specific internal energy instead of a linearly-spaced grid.
TabulatedBicubicFluidProperties
Fluid properties using bicubic interpolation on tabulated values provided
See TabulatedFluidProperties for more information on how to use this object.
warningwarning
Third order fluid property interpolations are not always monotonous. This is exacerbated by composition when using alternative variable sets such as (specific volume, specific internal energy). If this arises and is observed to remove the uniqueness of the numerical solution, lower order interpolation should be considered.
Input Parameters
- T_initial_guess400Temperature initial guess for Newton Method variable set conversion
Default:400
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature initial guess for Newton Method variable set conversion
- p_initial_guess200000Pressure initial guess for Newton Method variable set conversion
Default:200000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Pressure initial guess for Newton Method variable set conversion
- tolerance1e-08Tolerance for 2D Newton variable set conversion
Default:1e-08
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Tolerance for 2D Newton variable set conversion
Variable Set Conversions Newton Solve Parameters
- allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
- allow_imperfect_jacobiansFalsetrue to allow unimplemented property derivative terms to be set to zero for the AD API
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:true to allow unimplemented property derivative terms to be set to zero for the AD API
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Unit:(no unit assumed)
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
- force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
- force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
- force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup
- fp_typesingle-phase-fpType of the fluid property object
Default:single-phase-fp
C++ Type:FPType
Unit:(no unit assumed)
Controllable:No
Description:Type of the fluid property object
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- construct_pT_from_veFalseIf the lookup table (p, T) as functions of (v, e) should be constructed.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:If the lookup table (p, T) as functions of (v, e) should be constructed.
- construct_pT_from_vhFalseIf the lookup table (p, T) as functions of (v, h) should be constructed.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:If the lookup table (p, T) as functions of (v, h) should be constructed.
Variable Set Conversion Parameters
- e_maxMaximum specific internal energy for tabulated data.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum specific internal energy for tabulated data.
- e_minMinimum specific internal energy for tabulated data.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Minimum specific internal energy for tabulated data.
- out_of_bounds_behaviorthrowProperty evaluation behavior when evaluated outside the user-specified or tabulation-specified bounds
Default:throw
C++ Type:MooseEnum
Unit:(no unit assumed)
Options:ignore, throw, declare_invalid, set_to_closest_bound
Controllable:No
Description:Property evaluation behavior when evaluated outside the user-specified or tabulation-specified bounds
- pressure_max5e+07Maximum pressure for tabulated data.
Default:5e+07
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum pressure for tabulated data.
- pressure_min100000Minimum pressure for tabulated data.
Default:100000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Minimum pressure for tabulated data.
- temperature_max500Maximum temperature for tabulated data.
Default:500
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum temperature for tabulated data.
- temperature_min300Minimum temperature for tabulated data.
Default:300
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Minimum temperature for tabulated data.
- v_maxMaximum specific volume for tabulated data.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum specific volume for tabulated data.
- v_minMinimum specific volume for tabulated data.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Minimum specific volume for tabulated data.
Tabulation And Interpolation Bounds Parameters
- fluid_property_fileName of the csv file containing the tabulated fluid property data.
C++ Type:FileName
Unit:(no unit assumed)
Controllable:No
Description:Name of the csv file containing the tabulated fluid property data.
Tabulation File Read/Write Parameters
- num_T100Number of points to divide temperature range.
Default:100
C++ Type:unsigned int
Unit:(no unit assumed)
Controllable:No
Description:Number of points to divide temperature range.
- num_e100Number of points to divide specific internal energy range for (v,e) lookups.
Default:100
C++ Type:unsigned int
Unit:(no unit assumed)
Controllable:No
Description:Number of points to divide specific internal energy range for (v,e) lookups.
- num_p100Number of points to divide pressure range.
Default:100
C++ Type:unsigned int
Unit:(no unit assumed)
Controllable:No
Description:Number of points to divide pressure range.
- num_v100Number of points to divide specific volume range for (v,e) lookups.
Default:100
C++ Type:unsigned int
Unit:(no unit assumed)
Controllable:No
Description:Number of points to divide specific volume range for (v,e) lookups.
- use_log_grid_vFalseOption to use a base-10 logarithmically-spaced grid for specific volume instead of a linearly-spaced grid.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Option to use a base-10 logarithmically-spaced grid for specific volume instead of a linearly-spaced grid.
Tabulation And Interpolation Discretization Parameters
Input Files
- (modules/porous_flow/examples/reservoir_model/field_model.i)
- (modules/porous_flow/test/tests/sinks/injection_production_eg_outflowBC.i)
- (modules/fluid_properties/test/tests/tabulated/tabulated_v_e.i)
- (modules/porous_flow/test/tests/fluids/brine1_tabulated.i)
- (modules/porous_flow/examples/lava_lamp/2phase_convection.i)
- (modules/porous_flow/test/tests/fluidstate/water_vapor_tab.i)
- (modules/porous_flow/examples/tutorial/11_2D.i)
- (modules/porous_flow/examples/restart/gas_injection.i)
- (modules/porous_flow/test/tests/sinks/injection_production_eg.i)
- (modules/porous_flow/test/tests/fluidstate/theis_tabulated.i)
- (modules/porous_flow/examples/lava_lamp/1phase_convection.i)
- (modules/fluid_properties/test/tests/brine/brine_tabulated.i)
- (modules/porous_flow/examples/co2_intercomparison/1Dradial/1Dradial.i)
- (modules/porous_flow/examples/multiapp_fracture_flow/3dFracture/fracture_only_aperture_changing.i)
- (modules/fluid_properties/test/tests/tabulated/tabulated.i)
- (modules/porous_flow/examples/fluidflower/fluidflower.i)
- (modules/porous_flow/examples/reservoir_model/regular_grid.i)
- (modules/porous_flow/examples/restart/gas_injection_new_mesh.i)
(modules/porous_flow/examples/reservoir_model/field_model.i)
# Field model generated using geophysical modelling tool
[Mesh]
type = FileMesh
file = field.e
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 -9.81'
temperature_unit = Celsius
[]
[Problem]
# Variable porepressure has an initial condition despite the restart
allow_initial_conditions_with_restart = true
[]
[Variables]
[porepressure]
initial_condition = 20e6
[]
[]
[AuxVariables]
[temperature]
initial_condition = 50
[]
[xnacl]
initial_condition = 0.1
[]
[porosity]
family = MONOMIAL
order = CONSTANT
initial_from_file_var = poro
[]
[permx_md]
family = MONOMIAL
order = CONSTANT
initial_from_file_var = permX
[]
[permy_md]
family = MONOMIAL
order = CONSTANT
initial_from_file_var = permY
[]
[permz_md]
family = MONOMIAL
order = CONSTANT
initial_from_file_var = permZ
[]
[permx]
family = MONOMIAL
order = CONSTANT
[]
[permy]
family = MONOMIAL
order = CONSTANT
[]
[permz]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[permx]
type = ParsedAux
variable = permx
coupled_variables = permx_md
expression = '9.869233e-16*permx_md'
execute_on = initial
[]
[permy]
type = ParsedAux
variable = permy
coupled_variables = permy_md
expression = '9.869233e-16*permy_md'
execute_on = initial
[]
[permz]
type = ParsedAux
variable = permz
coupled_variables = permz_md
expression = '9.869233e-16*permz_md'
execute_on = initial
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
variable = porepressure
[]
[flux0]
type = PorousFlowFullySaturatedDarcyFlow
variable = porepressure
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = porepressure
number_fluid_phases = 1
number_fluid_components = 1
[]
[]
[FluidProperties]
[water]
type = Water97FluidProperties
[]
[watertab]
type = TabulatedBicubicFluidProperties
fp = water
save_file = false
[]
[]
[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
water_fp = watertab
[]
[porosity]
type = PorousFlowPorosityConst
porosity = porosity
[]
[permeability]
type = PorousFlowPermeabilityConstFromVar
perm_xx = permx
perm_yy = permy
perm_zz = permz
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1e2
end_time = 1e2
[]
[Outputs]
execute_on = 'initial timestep_end'
exodus = true
perf_graph = true
[]
(modules/porous_flow/test/tests/sinks/injection_production_eg_outflowBC.i)
# phase = 0 is liquid phase
# phase = 1 is gas phase
# fluid_component = 0 is water
# fluid_component = 1 is CO2
# Constant rates of water and CO2 injection into the left boundary
# 1D mesh
# The PorousFlowOutflowBCs remove the correct water and CO2 from the right boundary
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 20
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[AuxVariables]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[frac_water_in_liquid]
initial_condition = 1.0
[]
[frac_water_in_gas]
initial_condition = 0.0
[]
[water_kg_per_s]
[]
[co2_kg_per_s]
[]
[]
[AuxKernels]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = timestep_end
[]
[]
[Variables]
[pwater]
initial_condition = 20E6
[]
[pgas]
initial_condition = 21E6
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pwater
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pwater
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pgas
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = pgas
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas pwater'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[tabulated_water]
type = TabulatedBicubicFluidProperties
fp = true_water
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = water97_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated_11.csv
[]
[true_co2]
type = CO2FluidProperties
[]
[tabulated_co2]
type = TabulatedBicubicFluidProperties
fp = true_co2
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = co2_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = co2_tabulated_11.csv
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = 293.15
[]
[saturation_calculator]
type = PorousFlow2PhasePP
phase0_porepressure = pwater
phase1_porepressure = pgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'frac_water_in_liquid frac_water_in_gas'
[]
[water]
type = PorousFlowSingleComponentFluid
fp = tabulated_water
phase = 0
[]
[co2]
type = PorousFlowSingleComponentFluid
fp = tabulated_co2
phase = 1
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.2
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityBC
nw_phase = true
lambda = 2
s_res = 0.1
sum_s_res = 0.2
phase = 1
[]
[]
[BCs]
[water_injection]
type = PorousFlowSink
boundary = left
variable = pwater # pwater is associated with the water mass balance (fluid_component = 0 in its Kernels)
flux_function = -1E-5 # negative means a source, rather than a sink
[]
[co2_injection]
type = PorousFlowSink
boundary = left
variable = pgas # pgas is associated with the CO2 mass balance (fluid_component = 1 in its Kernels)
flux_function = -2E-5 # negative means a source, rather than a sink
[]
[right_water_component0]
type = PorousFlowOutflowBC
boundary = right
variable = pwater
mass_fraction_component = 0
save_in = water_kg_per_s
[]
[right_co2_component1]
type = PorousFlowOutflowBC
boundary = right
variable = pgas
mass_fraction_component = 1
save_in = co2_kg_per_s
[]
[]
[Preconditioning]
active = 'basic'
[basic]
type = SMP
full = true
petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu NONZERO 2'
[]
[preferred]
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
nl_abs_tol = 1E-10
nl_rel_tol = 1E-10
end_time = 1E5
[TimeStepper]
type = IterationAdaptiveDT
dt = 1E5
growth_factor = 1.1
[]
[]
[Postprocessors]
[water_kg_per_s]
type = NodalSum
boundary = right
variable = water_kg_per_s
[]
[co2_kg_per_s]
type = NodalSum
boundary = right
variable = co2_kg_per_s
[]
[]
[VectorPostprocessors]
[pps]
type = LineValueSampler
start_point = '0 0 0'
end_point = '20 0 0'
num_points = 20
sort_by = x
variable = 'pgas pwater saturation_gas'
[]
[]
[Outputs]
[out]
type = CSV
execute_on = final
[]
[]
(modules/fluid_properties/test/tests/tabulated/tabulated_v_e.i)
# Test thermophysical property calculations using TabulatedBiCubic/LinearFluidProperties.
# Calculations for density, internal energy and enthalpy using bicubic or bilinear
# interpolation of data generated using CO2FluidProperties.
[Mesh]
type = GeneratedMesh
dim = 2
# This test uses ElementalVariableValue postprocessors on specific
# elements, so element numbering needs to stay unchanged
allow_renumbering = false
[]
[AuxVariables]
[p]
family = MONOMIAL
order = CONSTANT
[]
[T]
family = MONOMIAL
order = CONSTANT
[]
[mu]
family = MONOMIAL
order = CONSTANT
[]
[s]
family = MONOMIAL
order = CONSTANT
[]
[cv]
family = MONOMIAL
order = CONSTANT
[]
[cp]
family = MONOMIAL
order = CONSTANT
[]
[c]
family = MONOMIAL
order = CONSTANT
[]
[k]
family = MONOMIAL
order = CONSTANT
[]
[g]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[pressure]
type = MaterialRealAux
variable = p
property = pressure
[]
[temperature]
type = MaterialRealAux
variable = T
property = temperature
[]
[viscosity]
type = MaterialRealAux
variable = mu
property = mu
[]
[s]
type = MaterialRealAux
variable = 's'
property = 's'
[]
[cv]
type = MaterialRealAux
variable = cv
property = cv
[]
[cp]
type = MaterialRealAux
variable = cp
property = cp
[]
[c]
type = MaterialRealAux
variable = c
property = c
[]
[thermal_conductivity]
type = MaterialRealAux
variable = k
property = k
[]
[g]
type = MaterialRealAux
variable = g
property = g
[]
[]
[FluidProperties]
[co2]
type = IdealGasFluidProperties
[]
[tabulated]
type = TabulatedBicubicFluidProperties
interpolated_properties = 'density enthalpy viscosity internal_energy k c cv cp entropy'
# Uncomment this to read the tabulation
# fluid_property_file = fluid_properties.csv
# Uncomment this to use the CO2 fluid properties above
# fp = 'co2'
# Uncomment this to write out a tabulation
# fluid_property_output_file = 'fluid_properties.csv'
# Enable the use of the (v,e) variables
construct_pT_from_ve = true
construct_pT_from_vh = true
out_of_bounds_behavior = 'set_to_closest_bound'
# Tabulation range
temperature_min = 280
temperature_max = 600
pressure_min = 1e5
pressure_max = 7e5
# Newton parameters
tolerance = 1e-8
T_initial_guess = 310
p_initial_guess = 1.8e5
[]
[]
[Materials]
[fp_mat_ve]
type = FluidPropertiesMaterialVE
v = 0.03108975251
e = -30797.6
fp = tabulated
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
[]
[Problem]
solve = false
[]
[Postprocessors]
[p]
type = ElementalVariableValue
elementid = 0
variable = p
[]
[T]
type = ElementalVariableValue
elementid = 0
variable = T
[]
[mu]
type = ElementalVariableValue
elementid = 0
variable = mu
[]
[s]
type = ElementalVariableValue
elementid = 0
variable = s
[]
[cv]
type = ElementalVariableValue
elementid = 0
variable = cv
[]
[cp]
type = ElementalVariableValue
elementid = 0
variable = cp
[]
[c]
type = ElementalVariableValue
elementid = 0
variable = c
[]
[k]
type = ElementalVariableValue
elementid = 0
variable = k
[]
[g]
type = ElementalVariableValue
elementid = 0
variable = g
[]
[]
[Outputs]
csv = true
file_base = tabulated_v_e_bilinear_out
execute_on = 'TIMESTEP_END'
[]
(modules/porous_flow/test/tests/fluids/brine1_tabulated.i)
# Test the density and viscosity calculated by the brine material using a
# TabulatedFluidProperties userobject for water
# Pressure 20 MPa
# Temperature 50C
# xnacl = 0.1047 (equivalent to 2.0 molality)
[Mesh]
type = GeneratedMesh
dim = 1
nx = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[]
[Variables]
[pp]
initial_condition = 20e6
[]
[]
[Kernels]
[dummy]
type = Diffusion
variable = pp
[]
[]
[AuxVariables]
[temp]
initial_condition = 50
[]
[xnacl]
initial_condition = 0.1047
[]
[]
[FluidProperties]
[water]
type = Water97FluidProperties
[]
[watertab]
type = TabulatedBicubicFluidProperties
fp = water
save_file = false
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = pp
[]
[brine]
type = PorousFlowBrine
water_fp = watertab
temperature_unit = Celsius
xnacl = 0.1047
phase = 0
[]
[]
[Executioner]
type = Steady
solve_type = Newton
[]
[Postprocessors]
[pressure]
type = ElementIntegralVariablePostprocessor
variable = pp
[]
[temperature]
type = ElementIntegralVariablePostprocessor
variable = temp
[]
[xnacl]
type = ElementIntegralVariablePostprocessor
variable = xnacl
[]
[density]
type = ElementIntegralMaterialProperty
mat_prop = 'PorousFlow_fluid_phase_density_qp0'
[]
[viscosity]
type = ElementIntegralMaterialProperty
mat_prop = 'PorousFlow_viscosity_qp0'
[]
[enthalpy]
type = ElementIntegralMaterialProperty
mat_prop = 'PorousFlow_fluid_phase_enthalpy_qp0'
[]
[internal_energy]
type = ElementIntegralMaterialProperty
mat_prop = 'PorousFlow_fluid_phase_internal_energy_qp0'
[]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = brine1
csv = true
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)
# Two phase density-driven convection of dissolved CO2 in brine
#
# Initially, the model has a gas phase at the top with a saturation of 0.29
# (which corresponds to an initial value of zi = 0.2).
# Diffusion of the dissolved CO2
# component from the saturated liquid to the unsaturated liquid below reduces the
# amount of CO2 in the gas phase. As the density of the CO2-saturated brine is greater
# than the unsaturated brine, a gravitational instability arises and density-driven
# convection of CO2-rich fingers descend into the unsaturated brine.
#
# The instability is seeded by a random perturbation to the porosity field.
# Mesh adaptivity is used to refine the mesh as the fingers form.
#
# Note: this model is computationally expensive, so should be run with multiple cores,
# preferably on a cluster.
[GlobalParams]
PorousFlowDictator = 'dictator'
gravity = '0 -9.81 0'
[]
[Adaptivity]
max_h_level = 2
marker = marker
initial_marker = initial
initial_steps = 2
[Indicators]
[indicator]
type = GradientJumpIndicator
variable = zi
[]
[]
[Markers]
[marker]
type = ErrorFractionMarker
indicator = indicator
refine = 0.8
[]
[initial]
type = BoxMarker
bottom_left = '0 1.95 0'
top_right = '2 2 0'
inside = REFINE
outside = DO_NOTHING
[]
[]
[]
[Mesh]
type = GeneratedMesh
dim = 2
ymax = 2
xmax = 2
ny = 40
nx = 40
bias_y = 0.95
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pgas
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pgas
[]
[diff0]
type = PorousFlowDispersiveFlux
fluid_component = 0
variable = pgas
disp_long = '0 0'
disp_trans = '0 0'
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = zi
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = zi
[]
[diff1]
type = PorousFlowDispersiveFlux
fluid_component = 1
variable = zi
disp_long = '0 0'
disp_trans = '0 0'
[]
[]
[AuxVariables]
[xnacl]
initial_condition = 0.01
[]
[saturation_gas]
order = FIRST
family = MONOMIAL
[]
[xco2l]
order = FIRST
family = MONOMIAL
[]
[density_liquid]
order = FIRST
family = MONOMIAL
[]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = 'timestep_end'
[]
[xco2l]
type = PorousFlowPropertyAux
variable = xco2l
property = mass_fraction
phase = 0
fluid_component = 1
execute_on = 'timestep_end'
[]
[density_liquid]
type = PorousFlowPropertyAux
variable = density_liquid
property = density
phase = 0
execute_on = 'timestep_end'
[]
[]
[Variables]
[pgas]
[]
[zi]
scaling = 1e4
[]
[]
[ICs]
[pressure]
type = FunctionIC
function = 10e6-9.81*1000*y
variable = pgas
[]
[zi]
type = BoundingBoxIC
variable = zi
x1 = 0
x2 = 2
y1 = 1.95
y2 = 2
inside = 0.2
outside = 0
[]
[porosity]
type = RandomIC
variable = porosity
min = 0.25
max = 0.275
seed = 0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
pc = 0
[]
[fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2
capillary_pressure = pc
[]
[]
[FluidProperties]
[co2sw]
type = CO2FluidProperties
[]
[co2]
type = TabulatedBicubicFluidProperties
fp = co2sw
[]
[brine]
type = BrineFluidProperties
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = '45'
[]
[brineco2]
type = PorousFlowFluidState
gas_porepressure = 'pgas'
z = 'zi'
temperature_unit = Celsius
xnacl = 'xnacl'
capillary_pressure = pc
fluid_state = fs
[]
[porosity]
type = PorousFlowPorosityConst
porosity = porosity
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-11 0 0 0 1e-11 0 0 0 1e-11'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
phase = 0
n = 2
s_res = 0.1
sum_s_res = 0.2
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
phase = 1
n = 2
s_res = 0.1
sum_s_res = 0.2
[]
[diffusivity]
type = PorousFlowDiffusivityConst
diffusion_coeff = '2e-9 2e-9 2e-9 2e-9'
tortuosity = '1 1'
[]
[]
[Preconditioning]
active = basic
[mumps_is_best_for_parallel_jobs]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[basic]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu NONZERO 2 '
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 1e6
nl_max_its = 25
l_max_its = 100
dtmax = 1e4
nl_abs_tol = 1e-6
[TimeStepper]
type = IterationAdaptiveDT
dt = 10
growth_factor = 2
cutback_factor = 0.5
[]
[]
[Functions]
[flux]
type = ParsedFunction
symbol_values = 'delta_xco2 dt'
symbol_names = 'dx dt'
expression = 'dx/dt'
[]
[]
[Postprocessors]
[total_co2_in_gas]
type = PorousFlowFluidMass
phase = 1
fluid_component = 1
[]
[total_co2_in_liquid]
type = PorousFlowFluidMass
phase = 0
fluid_component = 1
[]
[numdofs]
type = NumDOFs
[]
[delta_xco2]
type = ChangeOverTimePostprocessor
postprocessor = total_co2_in_liquid
[]
[dt]
type = TimestepSize
[]
[flux]
type = FunctionValuePostprocessor
function = flux
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
exodus = true
csv = true
[]
(modules/porous_flow/test/tests/fluidstate/water_vapor_tab.i)
# Tests correct calculation of properties in PorousFlowWaterVapor in the two-phase region
[Mesh]
type = GeneratedMesh
dim = 2
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pliq]
initial_condition = 1e6
[]
[h]
initial_condition = 8e5
scaling = 1e-3
[]
[]
[AuxVariables]
[pressure_gas]
order = CONSTANT
family = MONOMIAL
[]
[pressure_water]
order = CONSTANT
family = MONOMIAL
[]
[enthalpy_gas]
order = CONSTANT
family = MONOMIAL
[]
[enthalpy_water]
order = CONSTANT
family = MONOMIAL
[]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[saturation_water]
order = CONSTANT
family = MONOMIAL
[]
[density_water]
order = CONSTANT
family = MONOMIAL
[]
[density_gas]
order = CONSTANT
family = MONOMIAL
[]
[viscosity_water]
order = CONSTANT
family = MONOMIAL
[]
[viscosity_gas]
order = CONSTANT
family = MONOMIAL
[]
[temperature]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[enthalpy_water]
type = PorousFlowPropertyAux
variable = enthalpy_water
property = enthalpy
phase = 0
execute_on = 'initial timestep_end'
[]
[enthalpy_gas]
type = PorousFlowPropertyAux
variable = enthalpy_gas
property = enthalpy
phase = 1
execute_on = 'initial timestep_end'
[]
[pressure_water]
type = PorousFlowPropertyAux
variable = pressure_water
property = pressure
phase = 0
execute_on = 'initial timestep_end'
[]
[pressure_gas]
type = PorousFlowPropertyAux
variable = pressure_gas
property = pressure
phase = 1
execute_on = 'initial timestep_end'
[]
[saturation_water]
type = PorousFlowPropertyAux
variable = saturation_water
property = saturation
phase = 0
execute_on = 'initial timestep_end'
[]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = 'initial timestep_end'
[]
[density_water]
type = PorousFlowPropertyAux
variable = density_water
property = density
phase = 0
execute_on = 'initial timestep_end'
[]
[density_gas]
type = PorousFlowPropertyAux
variable = density_gas
property = density
phase = 1
execute_on = 'initial timestep_end'
[]
[viscosity_water]
type = PorousFlowPropertyAux
variable = viscosity_water
property = viscosity
phase = 0
execute_on = 'initial timestep_end'
[]
[viscosity_gas]
type = PorousFlowPropertyAux
variable = viscosity_gas
property = viscosity
phase = 1
execute_on = 'initial timestep_end'
[]
[temperature]
type = PorousFlowPropertyAux
variable = temperature
property = temperature
execute_on = 'initial timestep_end'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pliq h'
number_fluid_phases = 2
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureBC
pe = 1e5
lambda = 2
pc_max = 1e6
[]
[fs]
type = PorousFlowWaterVapor
water_fp = water
capillary_pressure = pc
[]
[]
[FluidProperties]
[water_true]
type = Water97FluidProperties
[]
[water]
type = TabulatedBicubicFluidProperties
fp = water_true
allow_fp_and_tabulation = true
fluid_property_file = fluid_properties_extended.csv
[]
[]
[Materials]
[watervapor]
type = PorousFlowFluidStateSingleComponent
porepressure = pliq
enthalpy = h
temperature_unit = Kelvin
capillary_pressure = pc
fluid_state = fs
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[relperm0]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[relperm1]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500
specific_heat_capacity = 1200
[]
[]
[Problem]
solve = false
[]
[Executioner]
type = Steady
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Postprocessors]
[density_water]
type = ElementAverageValue
variable = density_water
execute_on = 'initial timestep_end'
[]
[density_gas]
type = ElementAverageValue
variable = density_gas
execute_on = 'initial timestep_end'
[]
[viscosity_water]
type = ElementAverageValue
variable = viscosity_water
execute_on = 'initial timestep_end'
[]
[viscosity_gas]
type = ElementAverageValue
variable = viscosity_gas
execute_on = 'initial timestep_end'
[]
[enthalpy_water]
type = ElementAverageValue
variable = enthalpy_water
execute_on = 'initial timestep_end'
[]
[enthalpy_gas]
type = ElementAverageValue
variable = enthalpy_gas
execute_on = 'initial timestep_end'
[]
[sg]
type = ElementAverageValue
variable = saturation_gas
execute_on = 'initial timestep_end'
[]
[sw]
type = ElementAverageValue
variable = saturation_water
execute_on = 'initial timestep_end'
[]
[pwater]
type = ElementAverageValue
variable = pressure_water
execute_on = 'initial timestep_end'
[]
[pgas]
type = ElementAverageValue
variable = pressure_gas
execute_on = 'initial timestep_end'
[]
[temperature]
type = ElementAverageValue
variable = temperature
execute_on = 'initial timestep_end'
[]
[enthalpy]
type = ElementAverageValue
variable = h
execute_on = 'initial timestep_end'
[]
[liquid_mass]
type = PorousFlowFluidMass
phase = 0
execute_on = 'initial timestep_end'
[]
[vapor_mass]
type = PorousFlowFluidMass
phase = 1
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
file_base = water_vapor_twophase_tab
csv = true
execute_on = INITIAL
[]
(modules/porous_flow/examples/tutorial/11_2D.i)
# Two-phase borehole injection problem in RZ coordinates
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
xmin = 1.0
xmax = 10
bias_x = 1.4
ny = 3
ymin = -6
ymax = 6
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 -2 0'
top_right = '10 2 0'
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x<1.0001'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
coord_type = RZ
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pwater pgas T disp_r'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[]
[]
[GlobalParams]
displacements = 'disp_r disp_z'
gravity = '0 0 0'
biot_coefficient = 1.0
PorousFlowDictator = dictator
[]
[Variables]
[pwater]
initial_condition = 20E6
[]
[pgas]
initial_condition = 20.1E6
[]
[T]
initial_condition = 330
scaling = 1E-5
[]
[disp_r]
scaling = 1E-5
[]
[]
[Kernels]
[mass_water_dot]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pwater
[]
[flux_water]
type = PorousFlowAdvectiveFlux
fluid_component = 0
use_displaced_mesh = false
variable = pwater
[]
[vol_strain_rate_water]
type = PorousFlowMassVolumetricExpansion
fluid_component = 0
variable = pwater
[]
[mass_co2_dot]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pgas
[]
[flux_co2]
type = PorousFlowAdvectiveFlux
fluid_component = 1
use_displaced_mesh = false
variable = pgas
[]
[vol_strain_rate_co2]
type = PorousFlowMassVolumetricExpansion
fluid_component = 1
variable = pgas
[]
[energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = T
[]
[advection]
type = PorousFlowHeatAdvection
use_displaced_mesh = false
variable = T
[]
[conduction]
type = PorousFlowHeatConduction
use_displaced_mesh = false
variable = T
[]
[vol_strain_rate_heat]
type = PorousFlowHeatVolumetricExpansion
variable = T
[]
[grad_stress_r]
type = StressDivergenceRZTensors
temperature = T
variable = disp_r
eigenstrain_names = thermal_contribution
use_displaced_mesh = false
component = 0
[]
[poro_r]
type = PorousFlowEffectiveStressCoupling
variable = disp_r
use_displaced_mesh = false
component = 0
[]
[]
[AuxVariables]
[disp_z]
[]
[effective_fluid_pressure]
family = MONOMIAL
order = CONSTANT
[]
[mass_frac_phase0_species0]
initial_condition = 1 # all water in phase=0
[]
[mass_frac_phase1_species0]
initial_condition = 0 # no water in phase=1
[]
[sgas]
family = MONOMIAL
order = CONSTANT
[]
[swater]
family = MONOMIAL
order = CONSTANT
[]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_tt]
family = MONOMIAL
order = CONSTANT
[]
[stress_zz]
family = MONOMIAL
order = CONSTANT
[]
[porosity]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[effective_fluid_pressure]
type = ParsedAux
coupled_variables = 'pwater pgas swater sgas'
expression = 'pwater * swater + pgas * sgas'
variable = effective_fluid_pressure
[]
[swater]
type = PorousFlowPropertyAux
variable = swater
property = saturation
phase = 0
execute_on = timestep_end
[]
[sgas]
type = PorousFlowPropertyAux
variable = sgas
property = saturation
phase = 1
execute_on = timestep_end
[]
[stress_rr_aux]
type = RankTwoAux
variable = stress_rr
rank_two_tensor = stress
index_i = 0
index_j = 0
[]
[stress_tt]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_tt
index_i = 2
index_j = 2
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 1
index_j = 1
[]
[porosity]
type = PorousFlowPropertyAux
variable = porosity
property = porosity
execute_on = timestep_end
[]
[]
[BCs]
[pinned_top_bottom_r]
type = DirichletBC
variable = disp_r
value = 0
boundary = 'top bottom'
[]
[cavity_pressure_r]
type = Pressure
boundary = injection_area
variable = disp_r
postprocessor = constrained_effective_fluid_pressure_at_wellbore
use_displaced_mesh = false
[]
[cold_co2]
type = DirichletBC
boundary = injection_area
variable = T
value = 290 # injection temperature
use_displaced_mesh = false
[]
[constant_co2_injection]
type = PorousFlowSink
boundary = injection_area
variable = pgas
fluid_phase = 1
flux_function = -1E-4
use_displaced_mesh = false
[]
[outer_water_removal]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pwater
fluid_phase = 0
pt_vals = '0 1E9'
multipliers = '0 1E8'
PT_shift = 20E6
use_mobility = true
use_relperm = true
use_displaced_mesh = false
[]
[outer_co2_removal]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pgas
fluid_phase = 1
pt_vals = '0 1E9'
multipliers = '0 1E8'
PT_shift = 20.1E6
use_mobility = true
use_relperm = true
use_displaced_mesh = false
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[tabulated_water]
type = TabulatedBicubicFluidProperties
fp = true_water
temperature_min = 275
pressure_max = 1E8
fluid_property_output_file = water97_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated_11.csv
[]
[true_co2]
type = CO2FluidProperties
[]
[tabulated_co2]
type = TabulatedBicubicFluidProperties
fp = true_co2
temperature_min = 275
pressure_max = 1E8
fluid_property_output_file = co2_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = co2_tabulated_11.csv
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = T
[]
[saturation_calculator]
type = PorousFlow2PhasePP
phase0_porepressure = pwater
phase1_porepressure = pgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
[]
[water]
type = PorousFlowSingleComponentFluid
fp = tabulated_water
phase = 0
[]
[co2]
type = PorousFlowSingleComponentFluid
fp = tabulated_co2
phase = 1
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 4
s_res = 0.1
sum_s_res = 0.2
phase = 0
[]
[relperm_co2]
type = PorousFlowRelativePermeabilityBC
nw_phase = true
lambda = 2
s_res = 0.1
sum_s_res = 0.2
phase = 1
[]
[porosity]
type = PorousFlowPorosity
fluid = true
mechanical = true
thermal = true
porosity_zero = 0.1
reference_temperature = 330
reference_porepressure = 20E6
thermal_expansion_coeff = 15E-6 # volumetric
solid_bulk = 8E9 # unimportant since biot = 1
[]
[permeability_aquifer]
type = PorousFlowPermeabilityKozenyCarman
block = aquifer
poroperm_function = kozeny_carman_phi0
phi0 = 0.1
n = 2
m = 2
k0 = 1E-12
[]
[permeability_caps]
type = PorousFlowPermeabilityKozenyCarman
block = caps
poroperm_function = kozeny_carman_phi0
phi0 = 0.1
n = 2
m = 2
k0 = 1E-15
k_anisotropy = '1 0 0 0 1 0 0 0 0.1'
[]
[rock_thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2 0 0 0 2 0 0 0 2'
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 1100
density = 2300
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5E9
poissons_ratio = 0.0
[]
[strain]
type = ComputeAxisymmetricRZSmallStrain
eigenstrain_names = 'thermal_contribution initial_stress'
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = T
thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
eigenstrain_name = thermal_contribution
stress_free_temperature = 330
[]
[initial_strain]
type = ComputeEigenstrainFromInitialStress
initial_stress = '20E6 0 0 0 20E6 0 0 0 20E6'
eigenstrain_name = initial_stress
[]
[stress]
type = ComputeLinearElasticStress
[]
[effective_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[volumetric_strain]
type = PorousFlowVolumetricStrain
[]
[]
[Postprocessors]
[effective_fluid_pressure_at_wellbore]
type = PointValue
variable = effective_fluid_pressure
point = '1 0 0'
execute_on = timestep_begin
use_displaced_mesh = false
[]
[constrained_effective_fluid_pressure_at_wellbore]
type = FunctionValuePostprocessor
function = constrain_effective_fluid_pressure
execute_on = timestep_begin
[]
[]
[Functions]
[constrain_effective_fluid_pressure]
type = ParsedFunction
symbol_names = effective_fluid_pressure_at_wellbore
symbol_values = effective_fluid_pressure_at_wellbore
expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
[]
[]
[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 = 1E3
[TimeStepper]
type = IterationAdaptiveDT
dt = 1E3
growth_factor = 1.2
optimal_iterations = 10
[]
nl_abs_tol = 1E-7
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/examples/restart/gas_injection.i)
# 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]
uniform_refine = 1
[file]
type = FileMeshGenerator
file = gravityeq_out.e
[]
[translate]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0.1 0 0'
input = file
[]
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
symbol_values = injection_area
symbol_names = area
expression = '-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
[]
[]
[FluidProperties]
[brine]
type = BrineFluidProperties
[]
[methane]
type = MethaneFluidProperties
[]
[methane_tab]
type = TabulatedBicubicFluidProperties
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
checkpoint = true
[]
(modules/porous_flow/test/tests/sinks/injection_production_eg.i)
# phase = 0 is liquid phase
# phase = 1 is gas phase
# fluid_component = 0 is water
# fluid_component = 1 is CO2
# Constant rate of CO2 injection into the left boundary
# 1D mesh
# The PorousFlowPiecewiseLinearSinks remove the correct water and CO2 from the right boundary
# Note i take pretty big timesteps here so the system is quite nonlinear
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 20
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[AuxVariables]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[frac_water_in_liquid]
initial_condition = 1.0
[]
[frac_water_in_gas]
initial_condition = 0.0
[]
[]
[AuxKernels]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = timestep_end
[]
[]
[Variables]
[pwater]
initial_condition = 20E6
[]
[pgas]
initial_condition = 20.1E6
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pwater
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pwater
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pgas
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = pgas
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas pwater'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[tabulated_water]
type = TabulatedBicubicFluidProperties
fp = true_water
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = water97_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated_11.csv
[]
[true_co2]
type = CO2FluidProperties
[]
[tabulated_co2]
type = TabulatedBicubicFluidProperties
fp = true_co2
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = co2_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = co2_tabulated_11.csv
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = 293.15
[]
[saturation_calculator]
type = PorousFlow2PhasePP
phase0_porepressure = pwater
phase1_porepressure = pgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'frac_water_in_liquid frac_water_in_gas'
[]
[water]
type = PorousFlowSingleComponentFluid
fp = tabulated_water
phase = 0
[]
[co2]
type = PorousFlowSingleComponentFluid
fp = tabulated_co2
phase = 1
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.2
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityBC
nw_phase = true
lambda = 2
s_res = 0.1
sum_s_res = 0.2
phase = 1
[]
[]
[BCs]
[co2_injection]
type = PorousFlowSink
boundary = left
variable = pgas # pgas is associated with the CO2 mass balance (fluid_component = 1 in its Kernels)
flux_function = -1E-2 # negative means a source, rather than a sink
[]
[right_water]
type = PorousFlowPiecewiseLinearSink
boundary = right
# a sink of water, since the Kernels given to pwater are for fluid_component = 0 (the water)
variable = pwater
# this Sink is a function of liquid porepressure
# Also, all the mass_fraction, mobility and relperm are referenced to the liquid phase now
fluid_phase = 0
# Sink strength = (Pwater - 20E6)
pt_vals = '0 1E9'
multipliers = '0 1E9'
PT_shift = 20E6
# multiply Sink strength computed above by mass fraction of water at the boundary
mass_fraction_component = 0
# also multiply Sink strength by mobility of the liquid
use_mobility = true
# also multiply Sink strength by the relperm of the liquid
use_relperm = true
# also multiplly Sink strength by 1/L, where L is the distance to the fixed-porepressure external environment
flux_function = 10 # 1/L
[]
[right_co2]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pgas
fluid_phase = 1
pt_vals = '0 1E9'
multipliers = '0 1E9'
PT_shift = 20.1E6
mass_fraction_component = 1
use_mobility = true
use_relperm = true
flux_function = 10 # 1/L
[]
[]
[Preconditioning]
active = 'basic'
[basic]
type = SMP
full = true
petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu NONZERO 2'
[]
[preferred]
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
nl_abs_tol = 1E-13
nl_rel_tol = 1E-10
end_time = 1e4
[TimeStepper]
type = IterationAdaptiveDT
dt = 1E4
growth_factor = 1.1
[]
[]
[VectorPostprocessors]
[pps]
type = LineValueSampler
warn_discontinuous_face_values = false
start_point = '0 0 0'
end_point = '20 0 0'
num_points = 20
sort_by = x
variable = 'pgas pwater saturation_gas'
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/fluidstate/theis_tabulated.i)
# Two phase Theis problem: Flow from single source using WaterNCG fluidstate.
# Constant rate injection 2 kg/s
# 1D cylindrical mesh
# Initially, system has only a liquid phase, until enough gas is injected
# to form a gas phase, in which case the system becomes two phase.
# Note: this test is the same as theis.i, but uses the tabulated version of the CO2FluidProperties
[Mesh]
type = GeneratedMesh
dim = 1
nx = 80
xmax = 200
bias_x = 1.05
coord_type = RZ
rz_coord_axis = Y
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[AuxVariables]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[x1]
order = CONSTANT
family = MONOMIAL
[]
[y0]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = timestep_end
[]
[x1]
type = PorousFlowPropertyAux
variable = x1
property = mass_fraction
phase = 0
fluid_component = 1
execute_on = timestep_end
[]
[y0]
type = PorousFlowPropertyAux
variable = y0
property = mass_fraction
phase = 1
fluid_component = 0
execute_on = timestep_end
[]
[]
[Variables]
[pgas]
initial_condition = 20e6
[]
[zi]
initial_condition = 0
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pgas
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pgas
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = zi
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = zi
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
pc = 0
[]
[fs]
type = PorousFlowWaterNCG
water_fp = water
gas_fp = tabulated
capillary_pressure = pc
[]
[]
[FluidProperties]
[co2]
type = CO2FluidProperties
[]
[tabulated]
type = TabulatedBicubicFluidProperties
fp = co2
fluid_property_file = fluid_properties.csv
# We try to avoid using both, but some properties are not implemented in the tabulation
allow_fp_and_tabulation = true
# Test was design prior to bounds check
error_on_out_of_bounds = false
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = fluid_properties.csv
[]
[water]
type = Water97FluidProperties
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = 20
[]
[waterncg]
type = PorousFlowFluidState
gas_porepressure = pgas
z = zi
temperature_unit = Celsius
capillary_pressure = pc
fluid_state = fs
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.1
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 1
[]
[]
[BCs]
[rightwater]
type = DirichletBC
boundary = right
value = 20e6
variable = pgas
[]
[]
[DiracKernels]
[source]
type = PorousFlowSquarePulsePointSource
point = '0 0 0'
mass_flux = 2
variable = zi
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu NONZERO 2 1E-8 1E-10 20'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 8e2
[TimeStepper]
type = IterationAdaptiveDT
dt = 2
growth_factor = 2
[]
[]
[VectorPostprocessors]
[line]
type = LineValueSampler
warn_discontinuous_face_values = false
sort_by = x
start_point = '0 0 0'
end_point = '200 0 0'
num_points = 1000
variable = 'pgas zi x1 saturation_gas'
execute_on = 'timestep_end'
[]
[]
[Postprocessors]
[pgas]
type = PointValue
point = '1 0 0'
variable = pgas
[]
[sgas]
type = PointValue
point = '1 0 0'
variable = saturation_gas
[]
[zi]
type = PointValue
point = '1 0 0'
variable = zi
[]
[massgas]
type = PorousFlowFluidMass
fluid_component = 1
[]
[x1]
type = PointValue
point = '1 0 0'
variable = x1
[]
[y0]
type = PointValue
point = '1 0 0'
variable = y0
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
[csvout]
type = CSV
file_base = theis_tabulated_csvout
execute_on = timestep_end
execute_vector_postprocessors_on = final
[]
[]
(modules/porous_flow/examples/lava_lamp/1phase_convection.i)
# Two phase density-driven convection of dissolved CO2 in brine
#
# The model starts with CO2 in the liquid phase only. The CO2 diffuses into the brine.
# As the density of the CO2-saturated brine is greater
# than the unsaturated brine, a gravitational instability arises and density-driven
# convection of CO2-rich fingers descend into the unsaturated brine.
#
# The instability is seeded by a random perturbation to the porosity field.
# Mesh adaptivity is used to refine the mesh as the fingers form.
#
# Note: this model is computationally expensive, so should be run with multiple cores.
[GlobalParams]
PorousFlowDictator = 'dictator'
gravity = '0 -9.81 0'
[]
[Adaptivity]
max_h_level = 2
marker = marker
initial_marker = initial
initial_steps = 2
[Indicators]
[indicator]
type = GradientJumpIndicator
variable = zi
[]
[]
[Markers]
[marker]
type = ErrorFractionMarker
indicator = indicator
refine = 0.8
[]
[initial]
type = BoxMarker
bottom_left = '0 1.95 0'
top_right = '2 2 0'
inside = REFINE
outside = DO_NOTHING
[]
[]
[]
[Mesh]
type = GeneratedMesh
dim = 2
ymin = 1.5
ymax = 2
xmax = 2
ny = 20
nx = 40
bias_y = 0.95
[]
[AuxVariables]
[xnacl]
initial_condition = 0.01
[]
[saturation_gas]
order = FIRST
family = MONOMIAL
[]
[xco2l]
order = FIRST
family = MONOMIAL
[]
[density_liquid]
order = FIRST
family = MONOMIAL
[]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = 'timestep_end'
[]
[xco2l]
type = PorousFlowPropertyAux
variable = xco2l
property = mass_fraction
phase = 0
fluid_component = 1
execute_on = 'timestep_end'
[]
[density_liquid]
type = PorousFlowPropertyAux
variable = density_liquid
property = density
phase = 0
execute_on = 'timestep_end'
[]
[]
[Variables]
[pgas]
[]
[zi]
scaling = 1e4
[]
[]
[ICs]
[pressure]
type = FunctionIC
function = 10e6-9.81*1000*y
variable = pgas
[]
[zi]
type = ConstantIC
value = 0
variable = zi
[]
[porosity]
type = RandomIC
variable = porosity
min = 0.25
max = 0.275
seed = 0
[]
[]
[BCs]
[top]
type = DirichletBC
value = 0.04
variable = zi
boundary = top
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pgas
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pgas
[]
[diff0]
type = PorousFlowDispersiveFlux
fluid_component = 0
variable = pgas
disp_long = '0 0'
disp_trans = '0 0'
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = zi
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = zi
[]
[diff1]
type = PorousFlowDispersiveFlux
fluid_component = 1
variable = zi
disp_long = '0 0'
disp_trans = '0 0'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
pc = 0
[]
[fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2
capillary_pressure = pc
[]
[]
[FluidProperties]
[co2sw]
type = CO2FluidProperties
[]
[co2]
type = TabulatedBicubicFluidProperties
fp = co2sw
[]
[brine]
type = BrineFluidProperties
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = '45'
[]
[brineco2]
type = PorousFlowFluidState
gas_porepressure = 'pgas'
z = 'zi'
temperature_unit = Celsius
xnacl = 'xnacl'
capillary_pressure = pc
fluid_state = fs
[]
[porosity]
type = PorousFlowPorosityConst
porosity = porosity
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-11 0 0 0 1e-11 0 0 0 1e-11'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
phase = 0
n = 2
s_res = 0.1
sum_s_res = 0.2
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
phase = 1
n = 2
s_res = 0.1
sum_s_res = 0.2
[]
[diffusivity]
type = PorousFlowDiffusivityConst
diffusion_coeff = '2e-9 2e-9 2e-9 2e-9'
tortuosity = '1 1'
[]
[]
[Preconditioning]
active = basic
[mumps_is_best_for_parallel_jobs]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[basic]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu NONZERO 2 '
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 1e6
nl_max_its = 25
l_max_its = 100
dtmax = 1e4
nl_abs_tol = 1e-6
[TimeStepper]
type = IterationAdaptiveDT
dt = 100
growth_factor = 2
cutback_factor = 0.5
[]
[]
[Functions]
[flux]
type = ParsedFunction
symbol_values = 'delta_xco2 dt'
symbol_names = 'dx dt'
expression = 'dx/dt'
[]
[]
[Postprocessors]
[total_co2_in_gas]
type = PorousFlowFluidMass
phase = 1
fluid_component = 1
[]
[total_co2_in_liquid]
type = PorousFlowFluidMass
phase = 0
fluid_component = 1
[]
[numdofs]
type = NumDOFs
[]
[delta_xco2]
type = ChangeOverTimePostprocessor
postprocessor = total_co2_in_liquid
[]
[dt]
type = TimestepSize
[]
[flux]
type = FunctionValuePostprocessor
function = flux
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
exodus = true
csv = true
[]
(modules/fluid_properties/test/tests/brine/brine_tabulated.i)
# Test BrineFluidProperties calculations of density, viscosity and thermal
# conductivity with a TabulatedBiCubicFluidProperties water.
#
# Experimental density values from Pitzer et al, "Thermodynamic properties
# of aqueous sodium chloride solution", Journal of Physical and Chemical
# Reference Data, 13, 1-102 (1984)
#
# Experimental viscosity values from Phillips et al, "Viscosity of NaCl and
# other solutions up to 350C and 50MPa pressures", LBL-11586 (1980)
#
# Thermal conductivity values from Ozbek and Phillips, "Thermal conductivity of
# aqueous NaCl solutions from 20C to 330C", LBL-9086 (1980)
#
# --------------------------------------------------------------
# Pressure (Mpa) | 20 | 20 | 40
# Temperature (C) | 50 | 200 | 200
# NaCl molality (mol/kg) | 2 | 2 | 5
# NaCl mass fraction (kg/kg) | 0.1047 | 0.1047 | 0.2261
# --------------------------------------------------------------
# Expected values
# --------------------------------------------------------------
# Density (kg/m^3) | 1068.52 | 959.27 | 1065.58
# Viscosity (1e-6Pa.s) | 679.8 | 180.0 | 263.1
# Thermal conductivity (W/m/K) | 0.630 | 0.649 | 0.633
# --------------------------------------------------------------
# Calculated values
# --------------------------------------------------------------
# Density (kg/m^3) | 1067.18 | 958.68 | 1065.46
# Viscosity (1e-6 Pa.s) | 681.1 | 181.98 | 266.1
# Thermal conductivity (W/m/K) | 0.637 | 0.662 | 0.658
# --------------------------------------------------------------
#
# All results are within expected accuracy
[Mesh]
type = GeneratedMesh
dim = 2
nx = 3
ny = 1
xmax = 3
# This test uses ElementalVariableValue postprocessors on specific
# elements, so element numbering needs to stay unchanged
allow_renumbering = false
[]
[Variables]
[./dummy]
[../]
[]
[AuxVariables]
[./pressure]
family = MONOMIAL
order = CONSTANT
[../]
[./temperature]
family = MONOMIAL
order = CONSTANT
[../]
[./xnacl]
family = MONOMIAL
order = CONSTANT
[../]
[./density]
family = MONOMIAL
order = CONSTANT
[../]
[./enthalpy]
family = MONOMIAL
order = CONSTANT
[../]
[./internal_energy]
family = MONOMIAL
order = CONSTANT
[../]
[]
[Functions]
[./pic]
type = ParsedFunction
expression = 'if(x<2,20e6, 40e6)'
[../]
[./tic]
type = ParsedFunction
expression = 'if(x<1, 323.15, 473.15)'
[../]
[./xic]
type = ParsedFunction
expression = 'if(x<2,0.1047, 0.2261)'
[../]
[]
[ICs]
[./p_ic]
type = FunctionIC
function = pic
variable = pressure
[../]
[./t_ic]
type = FunctionIC
function = tic
variable = temperature
[../]
[./x_ic]
type = FunctionIC
function = xic
variable = xnacl
[../]
[]
[AuxKernels]
[./density]
type = MaterialRealAux
variable = density
property = density
[../]
[./enthalpy]
type = MaterialRealAux
variable = enthalpy
property = enthalpy
[../]
[./internal_energy]
type = MaterialRealAux
variable = internal_energy
property = e
[../]
[]
[FluidProperties]
[./water]
type = Water97FluidProperties
[../]
[./water_tab]
type = TabulatedBicubicFluidProperties
fp = water
save_file = false
[../]
[./brine]
type = BrineFluidProperties
water_fp = water_tab
[../]
[]
[Materials]
[./fp_mat]
type = MultiComponentFluidPropertiesMaterialPT
pressure = pressure
temperature = temperature
xmass = xnacl
fp = brine
[../]
[]
[Kernels]
[./diff]
type = Diffusion
variable = dummy
[../]
[]
[Executioner]
type = Steady
solve_type = NEWTON
[]
[Postprocessors]
[./density0]
type = ElementalVariableValue
variable = density
elementid = 0
[../]
[./density1]
type = ElementalVariableValue
variable = density
elementid = 1
[../]
[./density2]
type = ElementalVariableValue
variable = density
elementid = 2
[../]
[./enthalpy0]
type = ElementalVariableValue
variable = enthalpy
elementid = 0
[../]
[./enthalpy1]
type = ElementalVariableValue
variable = enthalpy
elementid = 1
[../]
[./enthalpy2]
type = ElementalVariableValue
variable = enthalpy
elementid = 2
[../]
[./e0]
type = ElementalVariableValue
variable = internal_energy
elementid = 0
[../]
[./e1]
type = ElementalVariableValue
variable = internal_energy
elementid = 1
[../]
[./e2]
type = ElementalVariableValue
variable = internal_energy
elementid = 2
[../]
[]
[Outputs]
csv = true
file_base = brine_out
[]
(modules/porous_flow/examples/co2_intercomparison/1Dradial/1Dradial.i)
# Intercomparison problem 3: Radial flow from an injection well
#
# From Pruess et al, Code intercomparison builds confidence in
# numerical simulation models for geologic disposal of CO2, Energy 29 (2004)
#
# A variation with zero salinity can be run by changing the initial condition
# of the AuxVariable xnacl
[Mesh]
type = GeneratedMesh
dim = 1
nx = 500
xmax = 10000
bias_x = 1.01
coord_type = 'RZ'
rz_coord_axis = Y
[]
[GlobalParams]
PorousFlowDictator = 'dictator'
gravity = '0 0 0'
[]
[AuxVariables]
[pressure_liquid]
order = CONSTANT
family = MONOMIAL
[]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[x1]
order = CONSTANT
family = MONOMIAL
[]
[y0]
order = CONSTANT
family = MONOMIAL
[]
[xnacl]
initial_condition = 0.15
[]
[]
[AuxKernels]
[pressure_liquid]
type = PorousFlowPropertyAux
variable = pressure_liquid
property = pressure
phase = 0
execute_on = 'timestep_end'
[]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = 'timestep_end'
[]
[x1]
type = PorousFlowPropertyAux
variable = x1
property = mass_fraction
phase = 0
fluid_component = 1
execute_on = 'timestep_end'
[]
[y0]
type = PorousFlowPropertyAux
variable = y0
property = mass_fraction
phase = 1
fluid_component = 0
execute_on = 'timestep_end'
[]
[]
[Variables]
[pgas]
initial_condition = 12e6
[]
[zi]
initial_condition = 0
scaling = 1e4
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pgas
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pgas
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = zi
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = zi
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi'
number_fluid_phases = 2
number_fluid_components = 3
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 5.099e-5
m = 0.457
sat_lr = 0.0
pc_max = 1e7
[]
[fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2
capillary_pressure = pc
[]
[]
[FluidProperties]
[co2sw]
type = CO2FluidProperties
[]
[co2]
type = TabulatedBicubicFluidProperties
fp = co2sw
[]
[water]
type = Water97FluidProperties
[]
[watertab]
type = TabulatedBicubicFluidProperties
fp = water
temperature_min = 273.15
temperature_max = 573.15
fluid_property_output_file = water_fluid_properties.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water_fluid_properties.csv
[]
[brine]
type = BrineFluidProperties
water_fp = watertab
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = '45'
[]
[brineco2]
type = PorousFlowFluidState
gas_porepressure = 'pgas'
z = 'zi'
temperature_unit = Celsius
xnacl = 'xnacl'
capillary_pressure = pc
fluid_state = fs
[]
[porosity]
type = PorousFlowPorosityConst
porosity = '0.12'
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityVG
m = 0.457
phase = 0
s_res = 0.3
sum_s_res = 0.35
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 1
s_res = 0.05
sum_s_res = 0.35
[]
[]
[BCs]
[rightwater]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
variable = pgas
use_mobility = true
PorousFlowDictator = dictator
fluid_phase = 0
multipliers = '0 1e9'
PT_shift = '12e6'
pt_vals = '0 1e9'
mass_fraction_component = 0
use_relperm = true
[]
[rightco2]
type = PorousFlowPiecewiseLinearSink
variable = zi
boundary = 'right'
use_mobility = true
PorousFlowDictator = dictator
fluid_phase = 1
multipliers = '0 1e9'
PT_shift = '12e6'
pt_vals = '0 1e9'
mass_fraction_component = 1
use_relperm = true
[]
[]
[DiracKernels]
[source]
type = PorousFlowSquarePulsePointSource
point = '0 0 0'
mass_flux = 1
variable = zi
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres bjacobi lu NONZERO'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 8.64e8
nl_max_its = 25
l_max_its = 100
dtmax = 5e6
[TimeStepper]
type = IterationAdaptiveDT
dt = 100
[]
[]
[VectorPostprocessors]
[vars]
type = NodalValueSampler
sort_by = x
variable = 'pgas zi xnacl'
execute_on = 'timestep_end'
outputs = spatial
[]
[auxvars]
type = ElementValueSampler
sort_by = x
variable = 'saturation_gas x1 y0'
execute_on = 'timestep_end'
outputs = spatial
[]
[]
[Postprocessors]
[pgas]
type = PointValue
point = '25.25 0 0'
variable = pgas
outputs = time
[]
[sgas]
type = PointValue
point = '25.25 0 0'
variable = saturation_gas
outputs = time
[]
[zi]
type = PointValue
point = '25.25 0 0'
variable = zi
outputs = time
[]
[massgas]
type = PorousFlowFluidMass
fluid_component = 1
outputs = time
[]
[x1]
type = PointValue
point = '25.25 0 0'
variable = x1
outputs = time
[]
[y0]
type = PointValue
point = '25.25 0 0'
variable = y0
outputs = time
[]
[xnacl]
type = PointValue
point = '25.25 0 0'
variable = xnacl
outputs = time
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
sync_times = '2.592e6 8.64e6 8.64e7 8.64e8'
[time]
type = CSV
[]
[spatial]
type = CSV
sync_only = true
[]
[]
(modules/porous_flow/examples/multiapp_fracture_flow/3dFracture/fracture_only_aperture_changing.i)
# Cold water injection into one side of the fracture network, and production from the other side
injection_rate = 10 # kg/s
[Mesh]
uniform_refine = 0
[cluster34]
type = FileMeshGenerator
file = 'Cluster_34.exo'
[]
[injection_node]
type = BoundingBoxNodeSetGenerator
input = cluster34
bottom_left = '-1000 0 -1000'
top_right = '1000 0.504 1000'
new_boundary = injection_node
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 -9.81E-6' # Note the value, because of pressure_unit
[]
[Variables]
[frac_P]
scaling = 1E6
[]
[frac_T]
initial_condition = 473
[]
[]
[ICs]
[frac_P]
type = FunctionIC
variable = frac_P
function = insitu_pp
[]
[]
[PorousFlowFullySaturated]
coupling_type = ThermoHydro
porepressure = frac_P
temperature = frac_T
fp = water
pressure_unit = MPa
[]
[Kernels]
[toMatrix]
type = PorousFlowHeatMassTransfer
variable = frac_T
v = transferred_matrix_T
transfer_coefficient = heat_transfer_coefficient
save_in = joules_per_s
[]
[]
[AuxVariables]
[heat_transfer_coefficient]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.0
[]
[transferred_matrix_T]
initial_condition = 473
[]
[joules_per_s]
[]
[normal_dirn_x]
family = MONOMIAL
order = CONSTANT
[]
[normal_dirn_y]
family = MONOMIAL
order = CONSTANT
[]
[normal_dirn_z]
family = MONOMIAL
order = CONSTANT
[]
[enclosing_element_normal_length]
family = MONOMIAL
order = CONSTANT
[]
[enclosing_element_normal_thermal_cond]
family = MONOMIAL
order = CONSTANT
[]
[aperture]
family = MONOMIAL
order = CONSTANT
[]
[perm_times_app]
family = MONOMIAL
order = CONSTANT
[]
[density]
family = MONOMIAL
order = CONSTANT
[]
[viscosity]
family = MONOMIAL
order = CONSTANT
[]
[insitu_pp]
[]
[]
[AuxKernels]
[normal_dirn_x_auxk]
type = PorousFlowElementNormal
variable = normal_dirn_x
component = x
[]
[normal_dirn_y]
type = PorousFlowElementNormal
variable = normal_dirn_y
component = y
[]
[normal_dirn_z]
type = PorousFlowElementNormal
variable = normal_dirn_z
component = z
[]
[heat_transfer_coefficient_auxk]
type = ParsedAux
variable = heat_transfer_coefficient
coupled_variables = 'enclosing_element_normal_length enclosing_element_normal_thermal_cond'
constant_names = h_s
constant_expressions = 1E3 # should be much bigger than thermal_conductivity / L ~ 1
expression = 'if(enclosing_element_normal_length = 0, 0, h_s * enclosing_element_normal_thermal_cond * 2 * enclosing_element_normal_length / (h_s * enclosing_element_normal_length * enclosing_element_normal_length + enclosing_element_normal_thermal_cond * 2 * enclosing_element_normal_length))'
[]
[aperture]
type = PorousFlowPropertyAux
variable = aperture
property = porosity
[]
[perm_times_app]
type = PorousFlowPropertyAux
variable = perm_times_app
property = permeability
row = 0
column = 0
[]
[density]
type = PorousFlowPropertyAux
variable = density
property = density
phase = 0
[]
[viscosity]
type = PorousFlowPropertyAux
variable = viscosity
property = viscosity
phase = 0
[]
[insitu_pp]
type = FunctionAux
execute_on = initial
variable = insitu_pp
function = insitu_pp
[]
[]
[BCs]
[inject_heat]
type = DirichletBC
boundary = injection_node
variable = frac_T
value = 373
[]
[]
[DiracKernels]
[inject_fluid]
type = PorousFlowPointSourceFromPostprocessor
mass_flux = ${injection_rate}
point = '58.8124 0.50384 74.7838'
variable = frac_P
[]
[withdraw_fluid]
type = PorousFlowPeacemanBorehole
SumQuantityUO = kg_out_uo
bottom_p_or_t = 10.6 # 1MPa + approx insitu at production point, to prevent aperture closing due to low porepressures
character = 1
line_length = 1
point_file = production.xyz
unit_weight = '0 0 0'
fluid_phase = 0
use_mobility = true
variable = frac_P
[]
[withdraw_heat]
type = PorousFlowPeacemanBorehole
SumQuantityUO = J_out_uo
bottom_p_or_t = 10.6 # 1MPa + approx insitu at production point, to prevent aperture closing due to low porepressures
character = 1
line_length = 1
point_file = production.xyz
unit_weight = '0 0 0'
fluid_phase = 0
use_mobility = true
use_enthalpy = true
variable = frac_T
[]
[]
[UserObjects]
[kg_out_uo]
type = PorousFlowSumQuantity
[]
[J_out_uo]
type = PorousFlowSumQuantity
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[water]
type = TabulatedBicubicFluidProperties
fp = true_water
temperature_min = 275 # K
temperature_max = 600
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = water97_tabulated.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated.csv
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosityLinear
porosity_ref = 1E-4 # fracture porosity = 1.0, but must include fracture aperture of 1E-4 at P = insitu_pp
P_ref = insitu_pp
P_coeff = 1E-3 # this is in metres/MPa, ie for P_ref = 1/P_coeff, the aperture becomes 1 metre
porosity_min = 1E-5
[]
[permeability]
type = PorousFlowPermeabilityKozenyCarman
k0 = 1E-15 # fracture perm = 1E-11 m^2, but must include fracture aperture of 1E-4
poroperm_function = kozeny_carman_phi0
m = 0
n = 3
phi0 = 1E-4
[]
[internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2700 # kg/m^3
specific_heat_capacity = 0 # basically no rock inside the fracture
[]
[aq_thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '0.6E-4 0 0 0 0.6E-4 0 0 0 0.6E-4' # thermal conductivity of water times fracture aperture. This should increase linearly with aperture, but is set constant in this model
[]
[]
[Functions]
[kg_rate]
type = ParsedFunction
symbol_values = 'dt kg_out'
symbol_names = 'dt kg_out'
expression = 'kg_out/dt'
[]
[insitu_pp]
type = ParsedFunction
expression = '10 - 0.847E-2 * z' # Approximate hydrostatic in MPa
[]
[]
[Postprocessors]
[dt]
type = TimestepSize
outputs = 'none'
[]
[kg_out]
type = PorousFlowPlotQuantity
uo = kg_out_uo
[]
[kg_per_s]
type = FunctionValuePostprocessor
function = kg_rate
[]
[J_out]
type = PorousFlowPlotQuantity
uo = J_out_uo
[]
[TK_out]
type = PointValue
variable = frac_T
point = '101.705 160.459 39.5722'
[]
[P_out]
type = PointValue
variable = frac_P
point = '101.705 160.459 39.5722'
[]
[P_in]
type = PointValue
variable = frac_P
point = '58.8124 0.50384 74.7838'
[]
[]
[VectorPostprocessors]
[heat_transfer_rate]
type = NodalValueSampler
outputs = none
sort_by = id
variable = joules_per_s
[]
[]
[Preconditioning]
[entire_jacobian]
type = SMP
full = true
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2 '
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
[TimeStepper]
type = IterationAdaptiveDT
dt = 1
optimal_iterations = 10
growth_factor = 1.5
[]
dtmax = 1E8
end_time = 1E8
nl_abs_tol = 1E-3
nl_max_its = 20
[]
[Outputs]
print_linear_residuals = false
csv = true
[ex]
type = Exodus
sync_times = '1 10 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 6100 6200 6300 6400 6500 6600 6700 6800 6900 7000 7100 7200 7300 7400 7500 7600 7700 7800 7900 8000 8100 8200 8300 8400 8500 8600 8700 8800 8900 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 30000 50000 70000 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000 2100000 2200000 2300000 2400000 2500000 2600000 2700000 2800000 2900000'
sync_only = true
[]
[]
(modules/fluid_properties/test/tests/tabulated/tabulated.i)
# Test thermophysical property calculations using TabulatedBiCubic/LinearFluidProperties.
# Calculations for density, internal energy and enthalpy using bicubic spline
# interpolation of data generated using CO2FluidProperties.
[Mesh]
type = GeneratedMesh
dim = 2
# This test uses ElementalVariableValue postprocessors on specific
# elements, so element numbering needs to stay unchanged
allow_renumbering = false
[]
[Variables]
[dummy]
[]
[]
[AuxVariables]
[pressure]
initial_condition = 2e6
family = MONOMIAL
order = CONSTANT
[]
[temperature]
initial_condition = 350
family = MONOMIAL
order = CONSTANT
[]
[rho]
family = MONOMIAL
order = CONSTANT
[]
[mu]
family = MONOMIAL
order = CONSTANT
[]
[e]
family = MONOMIAL
order = CONSTANT
[]
[h]
family = MONOMIAL
order = CONSTANT
[]
[s]
family = MONOMIAL
order = CONSTANT
[]
[cv]
family = MONOMIAL
order = CONSTANT
[]
[cp]
family = MONOMIAL
order = CONSTANT
[]
[c]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[rho]
type = MaterialRealAux
variable = rho
property = density
[]
[my]
type = MaterialRealAux
variable = mu
property = viscosity
[]
[internal_energy]
type = MaterialRealAux
variable = e
property = e
[]
[enthalpy]
type = MaterialRealAux
variable = h
property = h
[]
[entropy]
type = MaterialRealAux
variable = s
property = s
[]
[cv]
type = MaterialRealAux
variable = cv
property = cv
[]
[cp]
type = MaterialRealAux
variable = cp
property = cp
[]
[c]
type = MaterialRealAux
variable = c
property = c
[]
[]
[FluidProperties]
[co2]
type = CO2FluidProperties
[]
[tabulated]
type = TabulatedBicubicFluidProperties
fp = co2
interpolated_properties = 'density enthalpy viscosity internal_energy k c cv cp entropy'
# fluid_property_file = fluid_properties.csv
construct_pT_from_ve = false
construct_pT_from_vh = false
# Tabulation range
temperature_min = 280
temperature_max = 600
pressure_min = 1e5
pressure_max = 3e6
# Newton parameters
tolerance = 1e-8
T_initial_guess = 350
p_initial_guess = 1.5e5
[]
[]
[Materials]
[fp_mat]
type = FluidPropertiesMaterialPT
pressure = pressure
temperature = temperature
fp = tabulated
[]
[]
[Kernels]
[diff]
type = Diffusion
variable = dummy
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
[]
[Problem]
solve = false
[]
[Postprocessors]
[rho]
type = ElementalVariableValue
elementid = 0
variable = rho
[]
[mu]
type = ElementalVariableValue
elementid = 0
variable = mu
[]
[e]
type = ElementalVariableValue
elementid = 0
variable = e
[]
[h]
type = ElementalVariableValue
elementid = 0
variable = h
[]
[s]
type = ElementalVariableValue
elementid = 0
variable = s
[]
[cv]
type = ElementalVariableValue
elementid = 0
variable = cv
[]
[cp]
type = ElementalVariableValue
elementid = 0
variable = cp
[]
[c]
type = ElementalVariableValue
elementid = 0
variable = c
[]
[]
[Outputs]
csv = true
file_base = tabulated_out
execute_on = 'TIMESTEP_END'
perf_graph = true
[]
(modules/porous_flow/examples/fluidflower/fluidflower.i)
# FluidFlower International Benchmark study model
# CSIRO 2023
#
# This example can be used to reproduce the results presented by the
# CSIRO team as part of this benchmark study. See
# Green, C., Jackson, S.J., Gunning, J., Wilkins, A. and Ennis-King, J.,
# 2023. Modelling the FluidFlower: Insights from Characterisation and
# Numerical Predictions. Transport in Porous Media.
#
# This example takes a long time to run! The large density contrast
# between the gas phase CO2 and the water makes convergence very hard,
# so small timesteps must be taken during injection.
#
# This example uses a simplified mesh in order to be run during the
# automated testing. To reproduce the results of the benchmark study,
# replace the simple layered input mesh with the one located in the
# large_media submodule.
#
# The mesh file contains:
# - porosity as given by FluidFlower description
# - permeability as given by FluidFlower description
# - subdomain ids for each sand type
#
# The nominal thickness of the FluidFlower tank is 19mm. To keep masses consistent
# with the experiment, porosity and permeability are multiplied by the thickness
thickness = 0.019
#
# Properties associated with each sand type associated with mesh block ids
#
# block 0 - ESF (very fine sand)
sandESF = '0 10 20'
sandESF_pe = 1471.5
sandESF_krg = 0.09
sandESF_swi = 0.32
sandESF_krw = 0.71
sandESF_sgi = 0.14
# block 1 - C - Coarse lower
sandC = '1 21'
sandC_pe = 294.3
sandC_krg = 0.05
sandC_swi = 0.14
sandC_krw = 0.93
sandC_sgi = 0.1
# block 2 - D - Coarse upper
sandD = '2 22'
sandD_pe = 98.1
sandD_krg = 0.02
sandD_swi = 0.12
sandD_krw = 0.95
sandD_sgi = 0.08
# block 3 - E - Very Coarse lower
sandE = '3 13 23'
sandE_pe = 10
sandE_krg = 0.1
sandE_swi = 0.12
sandE_krw = 0.93
sandE_sgi = 0.06
# block 4 - F - Very Coarse upper
sandF = '4 14 24 34'
sandF_pe = 10
sandF_krg = 0.11
sandF_swi = 0.12
sandF_krw = 0.72
sandF_sgi = 0.13
# block 5 - G - Flush Zone
sandG = '5 15 35'
sandG_pe = 10
sandG_krg = 0.16
sandG_swi = 0.1
sandG_krw = 0.75
sandG_sgi = 0.06
# block 6 - Fault 1 - Heterogeneous
fault1 = '6 26'
fault1_pe = 10
fault1_krg = 0.16
fault1_swi = 0.1
fault1_krw = 0.75
fault1_sgi = 0.06
# block 7 - Fault 2 - Impermeable
# Note: this fault has been removed from the mesh (no elements in this region)
# block 8 - Fault 3 - Homogeneous
fault3 = '8'
fault3_pe = 10
fault3_krg = 0.16
fault3_swi = 0.1
fault3_krw = 0.75
fault3_sgi = 0.06
# Top layer
top_layer = '9'
# Boxes A, B an C used to report values (sg, sgr, xco2, etc)
boxA = '10 13 14 15 34 35'
boxB = '20 21 22 23 24 26'
boxC = '34 35'
# Furthermore, the seal sand unit in boxes A and B
seal_boxA = '10'
seal_boxB = '20'
# CO2 injection details:
# CO2 density ~1.8389 kg/m3 at 293.15 K, 1.01325e5 Pa
# Injection in Port (9, 3) for 5 hours.
# Injection in Port (17, 7) for 2:45 hours.
# Injection of 10 ml/min = 0.1666 ml/s = 1.666e-7 m3/s = ~3.06e-7 kg/s.
# Total mass of CO2 injected ~ 8.5g.
inj_rate = 3.06e-7
[Mesh]
[mesh]
type = FileMeshGenerator
file = 'fluidflower_test.e'
# file = '../../../../large_media/porous_flow/examples/fluidflower/fluidflower.e'
use_for_exodus_restart = true
[]
[]
[Debug]
show_var_residual_norms = true
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 -9.81 0'
temperature = temperature
log_extension = false
[]
[Variables]
[pgas]
family = MONOMIAL
order = CONSTANT
fv = true
[]
[z]
family = MONOMIAL
order = CONSTANT
fv = true
scaling = 1e4
[]
[]
[AuxVariables]
[xnacl]
family = MONOMIAL
order = CONSTANT
fv = true
initial_condition = 0.0055
[]
[temperature]
family = MONOMIAL
order = CONSTANT
fv = true
initial_condition = 20
[]
[porosity]
family = MONOMIAL
order = CONSTANT
fv = true
initial_from_file_var = porosity
[]
[porosity_times_thickness]
family = MONOMIAL
order = CONSTANT
fv = true
[]
[permeability]
family = MONOMIAL
order = CONSTANT
fv = true
initial_from_file_var = permeability
[]
[permeability_times_thickness]
family = MONOMIAL
order = CONSTANT
fv = true
[]
[saturation_water]
family = MONOMIAL
order = CONSTANT
[]
[saturation_gas]
family = MONOMIAL
order = CONSTANT
[]
[pressure_water]
family = MONOMIAL
order = CONSTANT
[]
[pc]
family = MONOMIAL
order = CONSTANT
[]
[x0_water]
order = CONSTANT
family = MONOMIAL
[]
[x0_gas]
order = CONSTANT
family = MONOMIAL
[]
[x1_water]
order = CONSTANT
family = MONOMIAL
[]
[x1_gas]
order = CONSTANT
family = MONOMIAL
[]
[density_water]
order = CONSTANT
family = MONOMIAL
[]
[density_gas]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[porosity_times_thickness]
type = ParsedAux
variable = porosity_times_thickness
coupled_variables = porosity
expression = 'porosity * ${thickness}'
execute_on = 'initial'
[]
[permeability_times_thickness]
type = ParsedAux
variable = permeability_times_thickness
coupled_variables = permeability
expression = 'permeability * ${thickness}'
execute_on = 'initial'
[]
[pressure_water]
type = ADPorousFlowPropertyAux
variable = pressure_water
property = pressure
phase = 0
execute_on = 'initial timestep_end'
[]
[saturation_water]
type = ADPorousFlowPropertyAux
variable = saturation_water
property = saturation
phase = 0
execute_on = 'initial timestep_end'
[]
[saturation_gas]
type = ADPorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = 'initial timestep_end'
[]
[density_water]
type = ADPorousFlowPropertyAux
variable = density_water
property = density
phase = 0
execute_on = 'initial timestep_end'
[]
[density_gas]
type = ADPorousFlowPropertyAux
variable = density_gas
property = density
phase = 1
execute_on = 'initial timestep_end'
[]
[x1_water]
type = ADPorousFlowPropertyAux
variable = x1_water
property = mass_fraction
phase = 0
fluid_component = 1
execute_on = 'initial timestep_end'
[]
[x1_gas]
type = ADPorousFlowPropertyAux
variable = x1_gas
property = mass_fraction
phase = 1
fluid_component = 1
execute_on = 'initial timestep_end'
[]
[x0_water]
type = ADPorousFlowPropertyAux
variable = x0_water
property = mass_fraction
phase = 0
fluid_component = 0
execute_on = 'initial timestep_end'
[]
[x0_gas]
type = ADPorousFlowPropertyAux
variable = x0_gas
property = mass_fraction
phase = 1
fluid_component = 0
execute_on = 'initial timestep_end'
[]
[pc]
type = ADPorousFlowPropertyAux
variable = pc
property = capillary_pressure
execute_on = 'initial timestep_end'
[]
[]
[FVKernels]
[mass0]
type = FVPorousFlowMassTimeDerivative
variable = pgas
fluid_component = 0
[]
[flux0]
type = FVPorousFlowAdvectiveFlux
variable = pgas
fluid_component = 0
[]
[diff0]
type = FVPorousFlowDispersiveFlux
variable = pgas
fluid_component = 0
disp_long = '0 0'
disp_trans = '0 0'
[]
[mass1]
type = FVPorousFlowMassTimeDerivative
variable = z
fluid_component = 1
[]
[flux1]
type = FVPorousFlowAdvectiveFlux
variable = z
fluid_component = 1
[]
[diff1]
type = FVPorousFlowDispersiveFlux
variable = z
fluid_component = 1
disp_long = '0 0'
disp_trans = '0 0'
[]
[]
[DiracKernels]
[injector1]
type = ConstantPointSource
point = '0.9 0.3 0'
value = ${inj_rate}
variable = z
[]
[injector2]
type = ConstantPointSource
point = '1.7 0.7 0'
value = ${inj_rate}
variable = z
[]
[]
[Controls]
[injection1]
type = ConditionalFunctionEnableControl
enable_objects = 'DiracKernels::injector1'
conditional_function = injection_schedule1
[]
[injection2]
type = ConditionalFunctionEnableControl
enable_objects = 'DiracKernels::injector2'
conditional_function = injection_schedule2
[]
[]
[Functions]
[initial_p]
type = ParsedFunction
symbol_names = 'p0 g H rho0'
symbol_values = '101.325e3 9.81 1.5 1002'
expression = 'p0 + rho0 * g * (H - y)'
[]
[injection_schedule1]
type = ParsedFunction
expression = 'if(t >= 0 & t <= 1.8e4, 1, 0)'
[]
[injection_schedule2]
type = ParsedFunction
expression = 'if(t >= 8.1e3 & t <= 1.8e4, 1, 0)'
[]
[]
[ICs]
[p]
type = FunctionIC
variable = pgas
function = initial_p
[]
[]
[FVBCs]
[pressure_top]
type = FVPorousFlowAdvectiveFluxBC
boundary = top
porepressure_value = 1.01325e5
variable = pgas
[]
[]
[FluidProperties]
[water]
type = Water97FluidProperties
[]
[watertab]
type = TabulatedBicubicFluidProperties
fp = water
save_file = false
pressure_min = 1e5
pressure_max = 1e6
temperature_min = 290
temperature_max = 300
num_p = 20
num_T = 10
[]
[co2]
type = CO2FluidProperties
[]
[co2tab]
type = TabulatedBicubicFluidProperties
fp = co2
save_file = false
pressure_min = 1e5
pressure_max = 1e6
temperature_min = 290
temperature_max = 300
num_p = 20
num_T = 10
[]
[brine]
type = BrineFluidProperties
water_fp = watertab
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas z'
number_fluid_phases = 2
number_fluid_components = 2
[]
[sandESF_pc]
type = PorousFlowCapillaryPressureBC
pe = ${sandESF_pe}
lambda = 2
block = ${sandESF}
pc_max = 1e4
sat_lr = ${sandESF_swi}
[]
[sandC_pc]
type = PorousFlowCapillaryPressureBC
pe = ${sandC_pe}
lambda = 2
block = ${sandC}
pc_max = 1e4
sat_lr = ${sandC_swi}
[]
[sandD_pc]
type = PorousFlowCapillaryPressureBC
pe = ${sandD_pe}
lambda = 2
block = ${sandD}
pc_max = 1e4
sat_lr = ${sandD_swi}
[]
[sandE_pc]
type = PorousFlowCapillaryPressureBC
pe = ${sandE_pe}
lambda = 2
block = ${sandE}
pc_max = 1e4
sat_lr = ${sandE_swi}
[]
[sandF_pc]
type = PorousFlowCapillaryPressureBC
pe = ${sandF_pe}
lambda = 2
block = ${sandF}
pc_max = 1e4
sat_lr = ${sandF_swi}
[]
[sandG_pc]
type = PorousFlowCapillaryPressureBC
pe = ${sandG_pe}
lambda = 2
block = ${sandG}
pc_max = 1e4
sat_lr = ${sandG_swi}
[]
[fault1_pc]
type = PorousFlowCapillaryPressureBC
pe = ${fault1_pe}
lambda = 2
block = ${fault1}
pc_max = 1e4
sat_lr = ${fault1_swi}
[]
[fault3_pc]
type = PorousFlowCapillaryPressureBC
pe = ${fault3_pe}
lambda = 2
block = ${fault3}
pc_max = 1e4
sat_lr = ${fault3_swi}
[]
[top_layer_pc]
type = PorousFlowCapillaryPressureConst
pc = 0
block = ${top_layer}
[]
[sandESF_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = sandESF_pc
[]
[sandC_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = sandC_pc
[]
[sandD_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = sandD_pc
[]
[sandE_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = sandE_pc
[]
[sandF_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = sandF_pc
[]
[sandG_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = sandG_pc
[]
[fault1_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = fault1_pc
[]
[fault3_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = fault3_pc
[]
[top_layer_fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2tab
capillary_pressure = top_layer_pc
[]
[]
[Materials]
[temperature]
type = ADPorousFlowTemperature
temperature = temperature
[]
[sandESF_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = sandESF_fs
capillary_pressure = sandESF_pc
block = ${sandESF}
[]
[sandC_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = sandC_fs
capillary_pressure = sandC_pc
block = ${sandC}
[]
[sandD_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = sandD_fs
capillary_pressure = sandD_pc
block = ${sandD}
[]
[sandE_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = sandE_fs
capillary_pressure = sandE_pc
block = ${sandE}
[]
[sandF_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = sandF_fs
capillary_pressure = sandF_pc
block = ${sandF}
[]
[sandG_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = sandG_fs
capillary_pressure = sandG_pc
block = ${sandG}
[]
[fault1_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = fault1_fs
capillary_pressure = fault1_pc
block = ${fault1}
[]
[fault3_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = fault3_fs
capillary_pressure = fault3_pc
block = ${fault3}
[]
[top_layer_brineco2]
type = ADPorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
fluid_state = top_layer_fs
capillary_pressure = top_layer_pc
block = ${top_layer}
[]
[porosity]
type = ADPorousFlowPorosityConst
porosity = porosity_times_thickness
[]
[permeability]
type = ADPorousFlowPermeabilityConstFromVar
perm_xx = permeability_times_thickness
perm_yy = permeability_times_thickness
perm_zz = permeability_times_thickness
[]
[diffcoeff]
type = ADPorousFlowDiffusivityConst
tortuosity = '1 1'
diffusion_coeff = '2e-9 2e-9 0 0'
[]
[sandESF_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${sandESF_swi}
sum_s_res = ${fparse sandESF_sgi + sandESF_swi}
scaling = ${sandESF_krw}
block = ${sandESF}
[]
[sandESF_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${sandESF_sgi}
sum_s_res = ${fparse sandESF_sgi + sandESF_swi}
scaling = ${sandESF_krg}
block = ${sandESF}
[]
[sandC_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${sandC_swi}
sum_s_res = ${fparse sandC_sgi + sandC_swi}
scaling = ${sandC_krw}
block = ${sandC}
[]
[sandC_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${sandC_sgi}
sum_s_res = ${fparse sandC_sgi + sandC_swi}
scaling = ${sandC_krg}
block = ${sandC}
[]
[sandD_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${sandD_swi}
sum_s_res = ${fparse sandD_sgi + sandD_swi}
scaling = ${sandD_krw}
block = ${sandD}
[]
[sandD_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${sandD_sgi}
sum_s_res = ${fparse sandD_sgi + sandD_swi}
scaling = ${sandD_krg}
block = ${sandD}
[]
[sandE_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${sandE_swi}
sum_s_res = ${fparse sandE_sgi + sandE_swi}
scaling = ${sandE_krw}
block = ${sandE}
[]
[sandE_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${sandE_sgi}
sum_s_res = ${fparse sandE_sgi + sandE_swi}
scaling = ${sandE_krg}
block = ${sandE}
[]
[sandF_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${sandF_swi}
sum_s_res = ${fparse sandF_sgi + sandF_swi}
scaling = ${sandF_krw}
block = ${sandF}
[]
[sandF_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${sandF_sgi}
sum_s_res = ${fparse sandF_sgi + sandF_swi}
scaling = ${sandF_krg}
block = ${sandF}
[]
[sandG_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${sandG_swi}
sum_s_res = ${fparse sandG_sgi + sandG_swi}
scaling = ${sandG_krw}
block = ${sandG}
[]
[sandG_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${sandG_sgi}
sum_s_res = ${fparse sandG_sgi + sandG_swi}
scaling = ${sandG_krg}
block = ${sandG}
[]
[fault1_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${fault1_swi}
sum_s_res = ${fparse fault1_sgi + fault1_swi}
scaling = ${fault1_krw}
block = ${fault1}
[]
[fault1_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${fault1_sgi}
sum_s_res = ${fparse fault1_sgi + fault1_swi}
scaling = ${fault1_krg}
block = ${fault1}
[]
[fault3_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
s_res = ${fault3_swi}
sum_s_res = ${fparse fault3_sgi + fault3_swi}
scaling = ${fault3_krw}
block = ${fault3}
[]
[fault3_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
s_res = ${fault3_sgi}
sum_s_res = ${fparse fault3_sgi + fault3_swi}
scaling = ${fault3_krg}
block = ${fault3}
[]
[top_layer_relperm0]
type = ADPorousFlowRelativePermeabilityBC
phase = 0
lambda = 2
block = ${top_layer}
[]
[top_layer_relperm1]
type = ADPorousFlowRelativePermeabilityBC
phase = 1
nw_phase = true
lambda = 2
block = ${top_layer}
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
petsc_options = '-ksp_snes_ew'
petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package -sub_pc_factor_shift_type'
petsc_options_value = 'gmres lu mumps NONZERO'
# petsc_options_iname = '-ksp_type -pc_type -pc_hypre_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels -ksp_gmres_restart'
# petsc_options_value = 'gmres hypre boomeramg lu NONZERO 4 301'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
dtmax = 60
start_time = 0
end_time = 4.32e5
nl_rel_tol = 1e-6
nl_abs_tol = 1e-8
nl_max_its = 15
l_tol = 1e-5
l_abs_tol = 1e-8
# line_search = none # Can be a useful option for this problem
[TimeSteppers]
[time]
type = FunctionDT
growth_factor = 2
cutback_factor_at_failure = 0.5
function = 'if(t<1.8e4, 2, if(t<3.6e4, 20, 60))'
[]
[]
[]
[Postprocessors]
[p_5_3]
type = PointValue
variable = pgas
point = '0.5 0.3 0'
execute_on = 'initial timestep_end'
[]
[p_5_3_w]
type = PointValue
variable = pressure_water
point = '0.5 0.3 0'
execute_on = 'initial timestep_end'
[]
[p_5_7]
type = PointValue
variable = pgas
point = '0.5 0.7 0'
execute_on = 'initial timestep_end'
[]
[p_5_7_w]
type = PointValue
variable = pressure_water
point = '0.5 0.7 0'
execute_on = 'initial timestep_end'
[]
[p_9_3]
type = PointValue
variable = pgas
point = '0.9 0.3 0'
execute_on = 'initial timestep_end'
[]
[p_9_3_w]
type = PointValue
variable = pressure_water
point = '0.9 0.3 0'
execute_on = 'initial timestep_end'
[]
[p_15_5]
type = PointValue
variable = pgas
point = '1.5 0.5 0'
execute_on = 'initial timestep_end'
[]
[p_15_5_w]
type = PointValue
variable = pressure_water
point = '1.5 0.5 0'
execute_on = 'initial timestep_end'
[]
[p_17_7]
type = PointValue
variable = pgas
point = '1.7 0.7 0'
execute_on = 'initial timestep_end'
[]
[p_17_7_w]
type = PointValue
variable = pressure_water
point = '1.7 0.7 0'
execute_on = 'initial timestep_end'
[]
[p_17_11]
type = PointValue
variable = pgas
point = '1.7 1.1 0'
execute_on = 'initial timestep_end'
[]
[p_17_11_w]
type = PointValue
variable = pressure_water
point = '1.7 1.1 0'
execute_on = 'initial timestep_end'
[]
[x0mass]
type = FVPorousFlowFluidMass
fluid_component = 0
phase = '0 1'
execute_on = 'initial timestep_end'
[]
[x1mass]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = '0 1'
execute_on = 'initial timestep_end'
[]
[x1gas]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = '1'
execute_on = 'initial timestep_end'
[]
[boxA]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = '0 1'
block = ${boxA}
execute_on = 'initial timestep_end'
[]
[imm_A_sandESF]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandESF_sgi}
block = 10
execute_on = 'initial timestep_end'
[]
[imm_A_sandE]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandE_sgi}
block = 13
execute_on = 'initial timestep_end'
[]
[imm_A_sandF]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandF_sgi}
block = '14 34'
execute_on = 'initial timestep_end'
[]
[imm_A_sandG]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandG_sgi}
block = '15 35'
execute_on = 'initial timestep_end'
[]
[imm_A]
type = LinearCombinationPostprocessor
pp_names = 'imm_A_sandESF imm_A_sandE imm_A_sandF imm_A_sandG'
pp_coefs = '1 1 1 1'
execute_on = 'initial timestep_end'
[]
[diss_A]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 0
block = ${boxA}
execute_on = 'initial timestep_end'
[]
[seal_A]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = '0 1'
block = ${seal_boxA}
execute_on = 'initial timestep_end'
[]
[boxB]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = '0 1'
block = ${boxB}
execute_on = 'initial timestep_end'
[]
[imm_B_sandESF]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandESF_sgi}
block = 20
execute_on = 'initial timestep_end'
[]
[imm_B_sandC]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandC_sgi}
block = 21
execute_on = 'initial timestep_end'
[]
[imm_B_sandD]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandD_sgi}
block = 22
execute_on = 'initial timestep_end'
[]
[imm_B_sandE]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandE_sgi}
block = 23
execute_on = 'initial timestep_end'
[]
[imm_B_sandF]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${sandF_sgi}
block = 24
execute_on = 'initial timestep_end'
[]
[imm_B_fault1]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 1
saturation_threshold = ${fault1_sgi}
block = 26
execute_on = 'initial timestep_end'
[]
[imm_B]
type = LinearCombinationPostprocessor
pp_names = 'imm_B_sandESF imm_B_sandC imm_B_sandD imm_B_sandE imm_B_sandF imm_B_fault1'
pp_coefs = '1 1 1 1 1 1'
execute_on = 'initial timestep_end'
[]
[diss_B]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = 0
block = ${boxB}
execute_on = 'initial timestep_end'
[]
[seal_B]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = '0 1'
block = ${seal_boxB}
execute_on = 'initial timestep_end'
[]
[boxC]
type = FVPorousFlowFluidMass
fluid_component = 1
phase = '0'
block = ${boxC}
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
# exodus = true
[csv]
type = CSV
[]
[]
(modules/porous_flow/examples/reservoir_model/regular_grid.i)
# SPE 10 comparative problem - model 1
# Data and description from https://www.spe.org/web/csp/datasets/set01.htm
# Simple input file that just establishes gravity equilibrium in the model
#
# Heterogeneous permeability is included by reading data from an external file
# using the PiecewiseMultilinear function, and saving that data to an elemental
# AuxVariable that is then used in PorousFlowPermeabilityConstFromVar
[Mesh]
type = GeneratedMesh
dim = 2
nx = 100
ny = 20
xmax = 762
ymax = 15.24
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 -9.81 0'
temperature_unit = Celsius
[]
[Variables]
[porepressure]
initial_condition = 20e6
[]
[]
[Functions]
[perm_md_fcn]
type = PiecewiseMultilinear
data_file = spe10_case1.data
[]
[]
[BCs]
[top]
type = DirichletBC
variable = porepressure
value = 20e6
boundary = top
[]
[]
[AuxVariables]
[temperature]
initial_condition = 50
[]
[xnacl]
initial_condition = 0.1
[]
[porosity]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.2
[]
[perm_md]
family = MONOMIAL
order = CONSTANT
[]
[perm]
family = MONOMIAL
order = CONSTANT
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
variable = porepressure
[]
[flux0]
type = PorousFlowFullySaturatedDarcyFlow
variable = porepressure
[]
[]
[AuxKernels]
[perm_md]
type = FunctionAux
function = perm_md_fcn
variable = perm_md
execute_on = initial
[]
[perm]
type = ParsedAux
variable = perm
coupled_variables = perm_md
expression = '9.869233e-16*perm_md'
execute_on = initial
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = porepressure
number_fluid_phases = 1
number_fluid_components = 1
[]
[]
[FluidProperties]
[water]
type = Water97FluidProperties
[]
[watertab]
type = TabulatedBicubicFluidProperties
fp = water
save_file = false
[]
[]
[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
water_fp = watertab
[]
[porosity]
type = PorousFlowPorosityConst
porosity = porosity
[]
[permeability]
type = PorousFlowPermeabilityConstFromVar
perm_xx = perm
perm_yy = perm
perm_zz = perm
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e5
nl_abs_tol = 1e-12
nl_rel_tol = 1e-06
steady_state_detection = true
steady_state_tolerance = 1e-12
[TimeStepper]
type = IterationAdaptiveDT
dt = 1e2
[]
[]
[Outputs]
execute_on = 'initial timestep_end'
exodus = true
perf_graph = true
[]
(modules/porous_flow/examples/restart/gas_injection_new_mesh.i)
# 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
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
symbol_values = injection_area
symbol_names = area
expression = '-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
[]
[]
[FluidProperties]
[brine]
type = BrineFluidProperties
[]
[methane]
type = MethaneFluidProperties
[]
[methane_tab]
type = TabulatedBicubicFluidProperties
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
[]