- uoPorousFlowSumQuantity user object name that holds the required information
C++ Type:UserObjectName
Description:PorousFlowSumQuantity user object name that holds the required information
PorousFlowPlotQuantity
This is used to record the total fluid (kg) or heat (J) flux that is produced by a PorousFlow DiracKernel in a time step. See polyline sinks for an extended discussion.
Input Parameters
- execute_onTIMESTEP_ENDThe list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Options:NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM, TRANSFER
Description:The list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM.
Optional 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
Options:
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).
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Options:
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Options:
Description:Set the enabled status of the MooseObject.
- force_preauxFalseForces the GeneralUserObject to be executed in PREAUX
Default:False
C++ Type:bool
Options:
Description:Forces the GeneralUserObject to be executed in PREAUX
- force_preicFalseForces the GeneralUserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Options:
Description:Forces the GeneralUserObject to be executed in PREIC during initial setup
- outputsVector of output names were you would like to restrict the output of variables(s) associated with this object
C++ Type:std::vector<OutputName>
Options:
Description:Vector of output names were you would like to restrict the output of variables(s) associated with this 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
Options:
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
Input Files
- (modules/porous_flow/test/tests/dirackernels/bh_except04.i)
- (modules/porous_flow/test/tests/dirackernels/pls03.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except11.i)
- (modules/porous_flow/test/tests/dirackernels/pls02.i)
- (modules/porous_flow/test/tests/dirackernels/bh05.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except01.i)
- (modules/porous_flow/examples/groundwater/ex02_abstraction.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except14.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except16.i)
- (modules/porous_flow/test/tests/dirackernels/bh02.i)
- (modules/porous_flow/test/tests/dirackernels/bh07.i)
- (modules/porous_flow/test/tests/dirackernels/bh03.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except12.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except10.i)
- (modules/porous_flow/test/tests/dirackernels/pls01.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except05.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except02.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except08.i)
- (modules/combined/examples/geochem-porous_flow/forge/porous_flow.i)
- (modules/combined/examples/geochem-porous_flow/geotes_2D/porous_flow.i)
- (modules/combined/examples/geochem-porous_flow/geotes_weber_tensleep/porous_flow.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except15.i)
- (modules/porous_flow/examples/groundwater/ex01.i)
- (modules/porous_flow/test/tests/dirackernels/pls03_action.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except03.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except07.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except13.i)
- (modules/porous_flow/examples/groundwater/ex02_steady_state.i)
- (modules/porous_flow/test/tests/dirackernels/bh04.i)
- (modules/porous_flow/test/tests/dirackernels/injection_production.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except09.i)
- (modules/porous_flow/test/tests/dirackernels/bh_except06.i)
(modules/porous_flow/test/tests/dirackernels/bh_except04.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
at_nodes = true # Needed to force exepected error
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
function_of = temperature
point_file = bh02.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/pls03.i)
# Test that the upwinding works correctly.
#
# A poly-line sink sits at the centre of the element.
# It has length=4 and weight=0.5, and extracts fluid
# at a constant rate of
# (1 * relative_permeability) kg.m^-1.s^-1
# Since it sits at the centre of the element, it extracts
# equally from each node, so the rate of extraction from
# each node is
# (0.5 * relative_permeability) kg.s^-1
# including the length and weight effects.
#
# There is no fluid flow.
#
# The initial conditions are such that all nodes have
# relative_permeability=0, except for one which has
# relative_permeaility = 1. Therefore, all nodes should
# remain at their initial porepressure, except the one.
#
# The porosity is 0.1, and the elemental volume is 2,
# so the fluid mass at the node in question = 0.2 * density / 4,
# where the 4 is the number of nodes in the element.
# In this simulation density = dens0 * exp(P / bulk), with
# dens0 = 100, and bulk = 20 MPa.
# The initial porepressure P0 = 10 MPa, so the final (after
# 1 second of simulation) is
# P(t=1) = 8.748592 MPa
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = 0
xmax = 2
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
#function = if((x<1)&(y<0.5),1E7,-1E7)
function = if((x<1)&(y>0.5),1E7,-1E7)
#function = if((x>1)&(y<0.5),1E7,-1E7)
#function = if((x>1)&(y>0.5),1E7,-1E7)
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pls_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e7
density0 = 100
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[relperm]
type = PorousFlowRelativePermeabilityFLAC
phase = 0
m = 2
s_res = 0.99
sum_s_res = 0.99
[]
[]
[DiracKernels]
[pls]
type = PorousFlowPolyLineSink
fluid_phase = 0
point_file = pls03.bh
use_relative_permeability = true
line_length = 4
SumQuantityUO = pls_total_outflow_mass
variable = pp
p_or_t_vals = '0 1E7'
fluxes = '1 1'
[]
[]
[Postprocessors]
[pls_report]
type = PorousFlowPlotQuantity
uo = pls_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p00]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[p01]
type = PointValue
variable = pp
point = '0 1 0'
execute_on = timestep_end
[]
[p20]
type = PointValue
variable = pp
point = '2 0 0'
execute_on = timestep_end
[]
[p21]
type = PointValue
variable = pp
point = '2 1 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 pls_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 1
dt = 1
solve_type = NEWTON
[]
[Outputs]
file_base = pls03
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except11.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
use_relative_permeability = true
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/pls02.i)
# fully-saturated situation with a poly-line sink with use_mobility=true
# The poly-line consists of 2 points, and has a length
# of 0.5. Each point is weighted with a weight of 0.1
# The PorousFlowPolyLineSink has
# p_or_t_vals = 0 1E7
# fluxes = 0 1
# so that for 0<=porepressure<=1E7
# base flux = porepressure * 1E-6 * mobility (measured in kg.m^-1.s^-1),
# and when multiplied by the poly-line length, and
# the weighting of each point, the mass flux is
# flux = porepressure * 0.5*E-8 * mobility (kg.s^-1).
#
# The fluid and matrix properties are:
# porosity = 0.1
# element volume = 8 m^3
# density = dens0 * exp(P / bulk), with bulk = 2E7
# initial porepressure P0 = 1E7
# viscosity = 0.2
# So, fluid mass = 0.8 * density (kg)
#
# The equation to solve is
# d(Mass)/dt = - porepressure * 0.5*E-8 * density / viscosity
#
# PorousFlow discretises time to conserve mass, so to march
# forward in time, we must solve
# Mass(dt) = Mass(0) - P * 0.5E-8 * density / viscosity * dt
# or
# 0.8 * dens0 * exp(P/bulk) = 0.8 * dens0 * exp(P0/bulk) - P * 0.5E-8 * density / viscosity * dt
# For the numbers written above this gives
# P(t=1) = 6.36947 MPa
# which is given precisely by MOOSE
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[pls_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e7
viscosity = 0.2
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[pls]
# This defines a sink that has strength
# f = L(P) * relperm * L_seg
# where
# L(P) is a piecewise-linear function of porepressure
# that is zero at pp=0 and 1 at pp=1E7
# relperm is the relative permeability of the fluid
# L_seg is the line-segment length associated with
# the Dirac points defined in the file pls02.bh
type = PorousFlowPolyLineSink
# Because the Variable for this Sink is pp, and pp is associated
# with the fluid-mass conservation equation, this sink is extracting
# fluid mass (and not heat energy or something else)
variable = pp
# The following specfies that the total fluid mass coming out of
# the porespace via this sink in this timestep should be recorded
# in the pls_total_outflow_mass UserObject
SumQuantityUO = pls_total_outflow_mass
# The following file defines the polyline geometry
# which is just two points in this particular example
point_file = pls02.bh
# Now define the piecewise-linear function, L
# First, we want L to be a function of porepressure (and not
# temperature or something else). The following means that
# p_or_t_vals should be intepreted by MOOSE as the zeroth-phase
# porepressure
function_of = pressure
fluid_phase = 0
# Second, define the piecewise-linear function, L
# The following means
# flux=0 when pp=0 (and also pp<0)
# flux=1 when pp=1E7 (and also pp>1E7)
# flux=linearly intepolated between pp=0 and pp=1E7
# When flux>0 this means a sink, while flux<0 means a source
p_or_t_vals = '0 1E7'
fluxes = '0 1'
# Finally, in this case we want to always multiply
# L by the fluid mobility (of the zeroth phase) and
# use that in the sink strength instead of the bare L
# computed above
use_mobility = true
[]
[]
[Postprocessors]
[pls_report]
type = PorousFlowPlotQuantity
uo = pls_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 pls_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 1
dt = 1
solve_type = NEWTON
[]
[Outputs]
file_base = pls02
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh05.i)
# unsaturated
# injection
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Functions]
[dts]
type = PiecewiseLinear
y = '500 500 1E1'
x = '4000 5000 6500'
[]
[]
[Variables]
[pp]
initial_condition = -2E5
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityFLAC
m = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
variable = pp
SumQuantityUO = borehole_total_outflow_mass
point_file = bh03.bh
fluid_phase = 0
bottom_p_or_t = 0
unit_weight = '0 0 0'
use_mobility = true
character = -1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 6500
solve_type = NEWTON
[TimeStepper]
type = FunctionDT
function = dts
[]
[]
[Outputs]
file_base = bh05
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except01.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 1
point_file = bh02.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/examples/groundwater/ex02_abstraction.i)
# Abstraction groundwater model. See groundwater_models.md for a detailed description
[Mesh]
[from_steady_state]
type = FileMeshGenerator
file = gold/ex02_steady_state_ex.e
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
function = steady_state_pp
[]
[]
[BCs]
[rainfall_recharge]
type = PorousFlowSink
boundary = zmax
variable = pp
flux_function = -1E-6 # recharge of 0.1mm/day = 1E-4m3/m2/day = 0.1kg/m2/day ~ 1E-6kg/m2/s
[]
[evapotranspiration]
type = PorousFlowHalfCubicSink
boundary = zmax
variable = pp
center = 0.0
cutoff = -5E4 # roots of depth 5m. 5m of water = 5E4 Pa
use_mobility = true
fluid_phase = 0
# Assume pan evaporation of 4mm/day = 4E-3m3/m2/day = 4kg/m2/day ~ 4E-5kg/m2/s
# Assume that if permeability was 1E-10m^2 and water table at topography then ET acts as pan strength
# Because use_mobility = true, then 4E-5 = maximum_flux = max * perm * density / visc = max * 1E-4, so max = 40
max = 40
[]
[]
[DiracKernels]
inactive = polyline_sink_borehole
[river]
type = PorousFlowPolyLineSink
SumQuantityUO = baseflow
point_file = ex02_river.bh
# Assume a perennial river.
# Assume the river has an incision depth of 1m and a stage height of 1.5m, and these are constant in time and uniform over the whole model. Hence, if groundwater head is 0.5m (5000Pa) there will be no baseflow and leakage.
p_or_t_vals = '-999995000 5000 1000005000'
# Assume the riverbed conductance, k_zz*density*river_segment_length*river_width/riverbed_thickness/viscosity = 1E-6*river_segment_length kg/Pa/s
fluxes = '-1E3 0 1E3'
variable = pp
[]
[horizontal_borehole]
type = PorousFlowPeacemanBorehole
SumQuantityUO = abstraction
bottom_p_or_t = -1E5
unit_weight = '0 0 -1E4'
character = 1.0
point_file = ex02.bh
variable = pp
[]
[polyline_sink_borehole]
type = PorousFlowPolyLineSink
SumQuantityUO = abstraction
fluxes = '-0.4 0 0.4'
p_or_t_vals = '-1E8 0 1E8'
point_file = ex02.bh
variable = pp
[]
[]
[Functions]
[steady_state_pp]
type = SolutionFunction
from_variable = pp
solution = steady_state_solution
[]
[baseflow_rate]
type = ParsedFunction
vars = 'baseflow_kg dt'
vals = 'baseflow_kg dt'
value = 'baseflow_kg / dt * 24.0 * 3600.0 / 400.0'
[]
[abstraction_rate]
type = ParsedFunction
vars = 'abstraction_kg dt'
vals = 'abstraction_kg dt'
value = 'abstraction_kg / dt * 24.0 * 3600.0'
[]
[]
[AuxVariables]
[ini_pp]
[]
[pp_change]
[]
[]
[AuxKernels]
[ini_pp]
type = FunctionAux
variable = ini_pp
function = steady_state_pp
execute_on = INITIAL
[]
[pp_change]
type = ParsedAux
variable = pp_change
args = 'pp ini_pp'
function = 'pp - ini_pp'
[]
[]
[PorousFlowUnsaturated]
fp = simple_fluid
porepressure = pp
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[]
[Materials]
[porosity_everywhere]
type = PorousFlowPorosityConst
porosity = 0.05
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst
block = 'top_aquifer bot_aquifer'
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-13'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst
block = aquitard
permeability = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
[]
[]
[UserObjects]
[steady_state_solution]
type = SolutionUserObject
execute_on = INITIAL
mesh = gold/ex02_steady_state_ex.e
timestep = LATEST
system_variables = pp
[]
[baseflow]
type = PorousFlowSumQuantity
[]
[abstraction]
type = PorousFlowSumQuantity
[]
[]
[Postprocessors]
[baseflow_kg]
type = PorousFlowPlotQuantity
uo = baseflow
outputs = 'none'
[]
[dt]
type = TimestepSize
outputs = 'none'
[]
[baseflow_l_per_m_per_day]
type = FunctionValuePostprocessor
function = baseflow_rate
[]
[abstraction_kg]
type = PorousFlowPlotQuantity
uo = abstraction
outputs = 'none'
[]
[abstraction_kg_per_day]
type = FunctionValuePostprocessor
function = abstraction_rate
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
# following 2 lines are not mandatory, but illustrate a popular preconditioner choice in groundwater models
petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap'
petsc_options_value = ' asm ilu 2 '
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 100
[TimeStepper]
type = FunctionDT
function = 'max(100, t)'
[]
end_time = 8.64E5 # 10 days
nl_abs_tol = 1E-11
[]
[Outputs]
print_linear_residuals = false
[ex]
type = Exodus
execute_on = final
[]
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/dirackernels/bh_except14.i)
# fully-saturated
# production
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02_huge.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except16.i)
# fully-saturated
# production
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
function_of = temperature
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh02.i)
# fully-saturated
# production
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
# Because the Variable for this Sink is pp, and pp is associated
# with the fluid-mass conservation equation, this sink is extracting
# fluid mass (and not heat energy or something else)
variable = pp
# The following specfies that the total fluid mass coming out of
# the porespace via this sink in this timestep should be recorded
# in the pls_total_outflow_mass UserObject
SumQuantityUO = borehole_total_outflow_mass
# The following file defines the polyline geometry
# which is just two points in this particular example
point_file = bh02.bh
# First, we want Peacemans f to be a function of porepressure (and not
# temperature or something else). So bottom_p_or_t is actually porepressure
function_of = pressure
fluid_phase = 0
# The bottomhole pressure
bottom_p_or_t = 0
# In this example there is no increase of the wellbore pressure
# due to gravity:
unit_weight = '0 0 0'
# PeacemanBoreholes should almost always have use_mobility = true
use_mobility = true
# This is a production wellbore (a sink of fluid that removes fluid from porespace)
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
[Outputs]
file_base = bh02
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh07.i)
# Comparison with analytical solution for cylindrically-symmetric situation
[Mesh]
type = FileMesh
file = bh07_input.e
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Functions]
[dts]
type = PiecewiseLinear
y = '1000 10000'
x = '100 1000'
[]
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[fflux]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pp
gravity = '0 0 0'
[]
[]
[BCs]
[fix_outer]
type = DirichletBC
boundary = perimeter
variable = pp
value = 1E7
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-11 0 0 0 1E-11 0 0 0 1E-11'
[]
[relperm]
type = PorousFlowRelativePermeabilityFLAC
m = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
variable = pp
SumQuantityUO = borehole_total_outflow_mass
point_file = bh07.bh
fluid_phase = 0
bottom_p_or_t = 0
unit_weight = '0 0 0'
use_mobility = true
re_constant = 0.1594 # use Chen and Zhang version
character = 2 # double the strength because bh07.bh only fills half the mesh
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
execute_on = 'initial timestep_end'
[]
[fluid_mass]
type = PorousFlowFluidMass
execute_on = 'initial timestep_end'
[]
[]
[VectorPostprocessors]
[pp]
type = LineValueSampler
variable = pp
start_point = '0 0 0'
end_point = '300 0 0'
sort_by = x
num_points = 300
execute_on = timestep_end
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 1E3
solve_type = NEWTON
[TimeStepper]
# get only marginally better results for smaller time steps
type = FunctionDT
function = dts
[]
[]
[Outputs]
file_base = bh07
[along_line]
type = CSV
execute_on = final
[]
[exodus]
type = Exodus
execute_on = 'initial final'
[]
[]
(modules/porous_flow/test/tests/dirackernels/bh03.i)
# fully-saturated
# injection
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 0
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
variable = pp
SumQuantityUO = borehole_total_outflow_mass
point_file = bh03.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 1E7
unit_weight = '0 0 0'
use_mobility = true
character = -1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
[Outputs]
file_base = bh03
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except12.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = does_not_exist
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except10.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
compute_internal_energy = false
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
use_mobility = true
use_internal_energy = true
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/pls01.i)
# fully-saturated situation with a poly-line sink at one
# of the nodes. Because there is no fluid flow, the
# other nodes should not experience any change in
# porepressure.
# The poly-line sink has length=2 and weight=0.1, and
# extracts fluid at a constant rate of 1 kg.m^-1.s^-1.
# Therefore, in 1 second it will have extracted a total
# of 0.2 kg.
# The porosity is 0.1, and the elemental volume is 2,
# so the fluid mass at the node in question = 0.2 * density / 4,
# where the 4 is the number of nodes in the element.
# In this simulation density = dens0 * exp(P / bulk), with
# dens0 = 100, and bulk = 20 MPa.
# The initial porepressure P0 = 10 MPa, so the final (after
# 1 second of simulation) is
# P(t=1) = 0.950879 MPa
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = 0
xmax = 2
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pls_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e7
density0 = 100
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[DiracKernels]
[pls]
type = PorousFlowPolyLineSink
fluid_phase = 0
point_file = pls01_21.bh
line_length = 2
SumQuantityUO = pls_total_outflow_mass
variable = pp
p_or_t_vals = '0 1E7'
fluxes = '1 1'
[]
[]
[Postprocessors]
[pls_report]
type = PorousFlowPlotQuantity
uo = pls_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p00]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[p01]
type = PointValue
variable = pp
point = '0 1 0'
execute_on = timestep_end
[]
[p20]
type = PointValue
variable = pp
point = '2 0 0'
execute_on = timestep_end
[]
[p21]
type = PointValue
variable = pp
point = '2 1 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 pls_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 1
dt = 1
solve_type = NEWTON
[]
[Outputs]
file_base = pls01
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except05.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
mass_fraction_component = 0
point_file = bh02.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except02.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
mass_fraction_component = 1
point_file = bh02.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except08.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
at_nodes = false # Needed to force expected error
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
use_mobility = true
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/combined/examples/geochem-porous_flow/forge/porous_flow.i)
# Input file modified from RobPodgorney version
# - 2D instead of 3D with different resolution. Effectively this means a 1m height of RobPodgorney aquifer is simulated. RobPodgorney total mass flux is 2.5kg/s meaning 0.25kg/s is appropriate here
# - Celsius instead of Kelvin
# - no use of PorousFlowPointEnthalpySourceFromPostprocessor since that is not yet merged into MOOSE: a DirichletBC is used instead
# - Use of PorousFlowFullySaturated instead of PorousFlowUnsaturated, and the save_component_rate_in feature to record the change in kg of each species at each node for passing to the Geochem simulation
# - MultiApps and Transfers to transfer information between this simulation and the aquifer_geochemistry.i simulation
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 225
ny = 200
xmin = -400
xmax = 500
ymin = -400
ymax = 400
[]
[injection_node]
input = gen
type = ExtraNodesetGenerator
new_boundary = injection_node
coord = '0 0 0'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[f_H]
initial_condition = 8.201229858451E-07
[]
[f_Na]
initial_condition = 2.281094143525E-03
[]
[f_K]
initial_condition = 2.305489507836E-04
[]
[f_Ca]
initial_condition = 5.818776782059E-04
[]
[f_Mg]
initial_condition = 1.539513498238E-07
[]
[f_SiO2]
initial_condition = 2.691822196469E-04
[]
[f_Al]
initial_condition = 4.457519474122E-08
[]
[f_Cl]
initial_condition = 4.744309776594E-03
[]
[f_SO4]
initial_condition = 9.516650880811E-06
[]
[f_HCO3]
initial_condition = 5.906126982324E-05
[]
[porepressure]
initial_condition = 20E6
[]
[temperature]
initial_condition = 220 # degC
scaling = 1E-6 # fluid enthalpy is roughly 1E6
[]
[]
[BCs]
[source_temperature]
type = DirichletBC
boundary = injection_node
variable = temperature
value = 70 # degC
[]
[]
[DiracKernels]
[inject_H]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 4.790385871045E-08
variable = f_H
[]
[inject_Na]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 7.586252963780E-07
variable = f_Na
[]
[inject_K]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 2.746517625125E-07
variable = f_K
[]
[inject_Ca]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 7.775129478597E-07
variable = f_Ca
[]
[inject_Mg]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 1.749872109005E-07
variable = f_Mg
[]
[inject_SiO2]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 4.100547515915E-06
variable = f_SiO2
[]
[inject_Al]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 2.502408592080E-08
variable = f_Al
[]
[inject_Cl]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 1.997260386272E-06
variable = f_Cl
[]
[inject_SO4]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 2.497372164191E-07
variable = f_SO4
[]
[inject_HCO3]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 5.003150992902E-06
variable = f_HCO3
[]
[inject_H2O]
type = PorousFlowPointSourceFromPostprocessor
point = ' 0 0 0'
mass_flux = 2.499865905987E-01
variable = porepressure
[]
[produce_H]
type = PorousFlowPeacemanBorehole
variable = f_H
SumQuantityUO = produced_mass_H
mass_fraction_component = 0
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_Na]
type = PorousFlowPeacemanBorehole
variable = f_Na
SumQuantityUO = produced_mass_Na
mass_fraction_component = 1
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_K]
type = PorousFlowPeacemanBorehole
variable = f_K
SumQuantityUO = produced_mass_K
mass_fraction_component = 2
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_Ca]
type = PorousFlowPeacemanBorehole
variable = f_Ca
SumQuantityUO = produced_mass_Ca
mass_fraction_component = 3
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_Mg]
type = PorousFlowPeacemanBorehole
variable = f_Mg
SumQuantityUO = produced_mass_Mg
mass_fraction_component = 4
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_SiO2]
type = PorousFlowPeacemanBorehole
variable = f_SiO2
SumQuantityUO = produced_mass_SiO2
mass_fraction_component = 5
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_Al]
type = PorousFlowPeacemanBorehole
variable = f_Al
SumQuantityUO = produced_mass_Al
mass_fraction_component = 6
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_Cl]
type = PorousFlowPeacemanBorehole
variable = f_Cl
SumQuantityUO = produced_mass_Cl
mass_fraction_component = 7
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_SO4]
type = PorousFlowPeacemanBorehole
variable = f_SO4
SumQuantityUO = produced_mass_SO4
mass_fraction_component = 8
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_HCO3]
type = PorousFlowPeacemanBorehole
variable = f_HCO3
SumQuantityUO = produced_mass_HCO3
mass_fraction_component = 9
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[produce_H2O]
type = PorousFlowPeacemanBorehole
variable = porepressure
SumQuantityUO = produced_mass_H2O
mass_fraction_component = 10
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[remove_heat_at_production_well]
type = PorousFlowPeacemanBorehole
variable = temperature
SumQuantityUO = produced_heat
point_file = production.bh
line_length = 1
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
use_enthalpy = true
character = 1
[]
[]
[UserObjects]
[produced_mass_H]
type = PorousFlowSumQuantity
[]
[produced_mass_Na]
type = PorousFlowSumQuantity
[]
[produced_mass_K]
type = PorousFlowSumQuantity
[]
[produced_mass_Ca]
type = PorousFlowSumQuantity
[]
[produced_mass_Mg]
type = PorousFlowSumQuantity
[]
[produced_mass_SiO2]
type = PorousFlowSumQuantity
[]
[produced_mass_Al]
type = PorousFlowSumQuantity
[]
[produced_mass_Cl]
type = PorousFlowSumQuantity
[]
[produced_mass_SO4]
type = PorousFlowSumQuantity
[]
[produced_mass_HCO3]
type = PorousFlowSumQuantity
[]
[produced_mass_H2O]
type = PorousFlowSumQuantity
[]
[produced_heat]
type = PorousFlowSumQuantity
[]
[]
[Postprocessors]
[heat_extracted]
type = PorousFlowPlotQuantity
uo = produced_heat
[]
[approx_production_temperature]
type = PointValue
point = '100 0 0'
variable = temperature
[]
[mass_extracted_H]
type = PorousFlowPlotQuantity
uo = produced_mass_H
execute_on = 'initial timestep_end'
[]
[mass_extracted_Na]
type = PorousFlowPlotQuantity
uo = produced_mass_Na
execute_on = 'initial timestep_end'
[]
[mass_extracted_K]
type = PorousFlowPlotQuantity
uo = produced_mass_K
execute_on = 'initial timestep_end'
[]
[mass_extracted_Ca]
type = PorousFlowPlotQuantity
uo = produced_mass_Ca
execute_on = 'initial timestep_end'
[]
[mass_extracted_Mg]
type = PorousFlowPlotQuantity
uo = produced_mass_Mg
execute_on = 'initial timestep_end'
[]
[mass_extracted_SiO2]
type = PorousFlowPlotQuantity
uo = produced_mass_SiO2
execute_on = 'initial timestep_end'
[]
[mass_extracted_Al]
type = PorousFlowPlotQuantity
uo = produced_mass_Al
execute_on = 'initial timestep_end'
[]
[mass_extracted_Cl]
type = PorousFlowPlotQuantity
uo = produced_mass_Cl
execute_on = 'initial timestep_end'
[]
[mass_extracted_SO4]
type = PorousFlowPlotQuantity
uo = produced_mass_SO4
execute_on = 'initial timestep_end'
[]
[mass_extracted_HCO3]
type = PorousFlowPlotQuantity
uo = produced_mass_HCO3
execute_on = 'initial timestep_end'
[]
[mass_extracted_H2O]
type = PorousFlowPlotQuantity
uo = produced_mass_H2O
execute_on = 'initial timestep_end'
[]
[mass_extracted]
type = LinearCombinationPostprocessor
pp_names = 'mass_extracted_H mass_extracted_Na mass_extracted_K mass_extracted_Ca mass_extracted_Mg mass_extracted_SiO2 mass_extracted_Al mass_extracted_Cl mass_extracted_SO4 mass_extracted_HCO3 mass_extracted_H2O'
pp_coefs = '1 1 1 1 1 1 1 1 1 1 1'
execute_on = 'initial timestep_end'
[]
[dt]
type = TimestepSize
execute_on = 'timestep_begin'
[]
[]
[Modules]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 2E-4
bulk_modulus = 2E9
viscosity = 1E-3
density0 = 980
cv = 4000.0
cp = 4000.0
porepressure_coefficient = 0
[]
[]
[]
[PorousFlowFullySaturated]
coupling_type = ThermoHydro
porepressure = porepressure
temperature = temperature
mass_fraction_vars = 'f_H f_Na f_K f_Ca f_Mg f_SiO2 f_Al f_Cl f_SO4 f_HCO3'
save_component_rate_in = 'rate_H rate_Na rate_K rate_Ca rate_Mg rate_SiO2 rate_Al rate_Cl rate_SO4 rate_HCO3 rate_H2O' # change in kg at every node / dt
fp = the_simple_fluid
temperature_unit = Celsius
[]
[AuxVariables]
[rate_H]
[]
[rate_Na]
[]
[rate_K]
[]
[rate_Ca]
[]
[rate_Mg]
[]
[rate_SiO2]
[]
[rate_Al]
[]
[rate_Cl]
[]
[rate_SO4]
[]
[rate_HCO3]
[]
[rate_H2O]
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.01
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.5 0 0 0 2.5 0 0 0 2.5'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
density = 2750.0
specific_heat_capacity = 900.0
[]
[relperm] # this is necessary because of the use_mobility = true used in the PorousFlowPeacemanBorehole objects
type = PorousFlowRelativePermeabilityConst
at_nodes = true
phase = 0
[]
[]
[Preconditioning]
active = typically_efficient
[typically_efficient]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = ' hypre boomeramg'
[]
[strong]
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 ilu NONZERO 2'
[]
[probably_too_strong]
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 = 31536000 #1 year
[TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 500
[]
[]
[Outputs]
exodus = true
csv = true
[]
[MultiApps]
[react]
type = TransientMultiApp
input_files = aquifer_geochemistry.i
clone_master_mesh = true
execute_on = 'timestep_end'
[]
[]
[Transfers]
[changes_due_to_flow]
type = MultiAppCopyTransfer
direction = to_multiapp
source_variable = 'rate_H rate_Na rate_K rate_Ca rate_Mg rate_SiO2 rate_Al rate_Cl rate_SO4 rate_HCO3 rate_H2O temperature'
variable = 'pf_rate_H pf_rate_Na pf_rate_K pf_rate_Ca pf_rate_Mg pf_rate_SiO2 pf_rate_Al pf_rate_Cl pf_rate_SO4 pf_rate_HCO3 pf_rate_H2O temperature'
multi_app = react
[]
[massfrac_from_geochem]
type = MultiAppCopyTransfer
direction = from_multiapp
source_variable = 'massfrac_H massfrac_Na massfrac_K massfrac_Ca massfrac_Mg massfrac_SiO2 massfrac_Al massfrac_Cl massfrac_SO4 massfrac_HCO3'
variable = 'f_H f_Na f_K f_Ca f_Mg f_SiO2 f_Al f_Cl f_SO4 f_HCO3'
multi_app = react
[]
[]
(modules/combined/examples/geochem-porous_flow/geotes_2D/porous_flow.i)
# PorousFlow simulation of injection and production in a 2D aquifer
# Much of this file is standard porous-flow stuff. The unusual aspects are:
# - transfer of the rates of changes of each species (kg/s) to the aquifer_geochemistry.i simulation. This is achieved by saving these changes from the PorousFlowMassTimeDerivative residuals
# - transfer of the temperature field to the aquifer_geochemistry.i simulation
# Interesting behaviour can be simulated by this file without its "parent" simulation, exchanger.i. exchanger.i provides mass-fractions injected via the injection_rate_massfrac_* variables, but since these are more-or-less constant throughout the duration of the exchanger.i simulation, the initial_conditions specified below may be used. Similar, exchanger.i provides injection_temperature, but that is also constant.
injection_rate = -1.0 # kg/s/m, negative because injection as a source
production_rate = 1.0 # kg/s/m
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 14 # for better resolution, use 56 or 112
ny = 8 # for better resolution, use 32 or 64
xmin = -70
xmax = 70
ymin = -40
ymax = 40
[]
[injection_node]
input = gen
type = ExtraNodesetGenerator
new_boundary = injection_node
coord = '-30 0 0'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[f0]
initial_condition = 0.002285946
[]
[f1]
initial_condition = 0.0035252
[]
[f2]
initial_condition = 1.3741E-05
[]
[porepressure]
initial_condition = 2E6
[]
[temperature]
initial_condition = 50
scaling = 1E-6 # fluid enthalpy is roughly 1E6
[]
[]
[BCs]
[injection_temperature]
type = MatchedValueBC
variable = temperature
v = injection_temperature
boundary = injection_node
[]
[]
[DiracKernels]
[inject_Na]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
line_length = 1.0
multiplying_var = injection_rate_massfrac_Na
point_file = injection.bh
variable = f0
[]
[inject_Cl]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
line_length = 1.0
multiplying_var = injection_rate_massfrac_Cl
point_file = injection.bh
variable = f1
[]
[inject_SiO2]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
line_length = 1.0
multiplying_var = injection_rate_massfrac_SiO2
point_file = injection.bh
variable = f2
[]
[inject_H2O]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
line_length = 1.0
multiplying_var = injection_rate_massfrac_H2O
point_file = injection.bh
variable = porepressure
[]
[produce_Na]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Na
fluxes = ${production_rate}
p_or_t_vals = 0.0
line_length = 1.0
mass_fraction_component = 0
point_file = production.bh
variable = f0
[]
[produce_Cl]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Cl
fluxes = ${production_rate}
p_or_t_vals = 0.0
line_length = 1.0
mass_fraction_component = 1
point_file = production.bh
variable = f1
[]
[produce_SiO2]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_SiO2
fluxes = ${production_rate}
p_or_t_vals = 0.0
line_length = 1.0
mass_fraction_component = 2
point_file = production.bh
variable = f2
[]
[produce_H2O]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_H2O
fluxes = ${production_rate}
p_or_t_vals = 0.0
line_length = 1.0
mass_fraction_component = 3
point_file = production.bh
variable = porepressure
[]
[produce_heat]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_heat
fluxes = ${production_rate}
p_or_t_vals = 0.0
line_length = 1.0
use_enthalpy = true
point_file = production.bh
variable = temperature
[]
[]
[UserObjects]
[injected_mass]
type = PorousFlowSumQuantity
[]
[produced_mass_Na]
type = PorousFlowSumQuantity
[]
[produced_mass_Cl]
type = PorousFlowSumQuantity
[]
[produced_mass_SiO2]
type = PorousFlowSumQuantity
[]
[produced_mass_H2O]
type = PorousFlowSumQuantity
[]
[produced_heat]
type = PorousFlowSumQuantity
[]
[]
[Postprocessors]
[dt]
type = TimestepSize
execute_on = TIMESTEP_BEGIN
[]
[tot_kg_injected_this_timestep]
type = PorousFlowPlotQuantity
uo = injected_mass
[]
[kg_Na_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Na
[]
[kg_Cl_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Cl
[]
[kg_SiO2_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_SiO2
[]
[kg_H2O_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_H2O
[]
[mole_rate_Na_produced]
type = FunctionValuePostprocessor
function = moles_Na
[]
[mole_rate_Cl_produced]
type = FunctionValuePostprocessor
function = moles_Cl
[]
[mole_rate_SiO2_produced]
type = FunctionValuePostprocessor
function = moles_SiO2
[]
[mole_rate_H2O_produced]
type = FunctionValuePostprocessor
function = moles_H2O
[]
[heat_joules_extracted_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_heat
[]
[production_temperature]
type = PointValue
point = '30 0 0'
variable = temperature
[]
[]
[Functions]
[moles_Na]
type = ParsedFunction
vars = 'kg_Na dt'
vals = 'kg_Na_produced_this_timestep dt'
value = 'kg_Na * 1000 / 22.9898 / dt'
[]
[moles_Cl]
type = ParsedFunction
vars = 'kg_Cl dt'
vals = 'kg_Cl_produced_this_timestep dt'
value = 'kg_Cl * 1000 / 35.453 / dt'
[]
[moles_SiO2]
type = ParsedFunction
vars = 'kg_SiO2 dt'
vals = 'kg_SiO2_produced_this_timestep dt'
value = 'kg_SiO2 * 1000 / 60.0843 / dt'
[]
[moles_H2O]
type = ParsedFunction
vars = 'kg_H2O dt'
vals = 'kg_H2O_produced_this_timestep dt'
value = 'kg_H2O * 1000 / 18.0152 / dt'
[]
[]
[Modules]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0
bulk_modulus = 2E9
viscosity = 1E-3
density0 = 1000
cv = 4000.0
cp = 4000.0
[]
[]
[]
[PorousFlowFullySaturated]
coupling_type = ThermoHydro
porepressure = porepressure
temperature = temperature
mass_fraction_vars = 'f0 f1 f2'
save_component_rate_in = 'rate_Na rate_Cl rate_SiO2 rate_H2O' # change in kg at every node / dt
fp = the_simple_fluid
temperature_unit = Celsius
[]
[AuxVariables]
[injection_temperature]
initial_condition = 200
[]
[injection_rate_massfrac_Na]
initial_condition = 0.002285946
[]
[injection_rate_massfrac_Cl]
initial_condition = 0.0035252
[]
[injection_rate_massfrac_SiO2]
initial_condition = 1.3741E-05
[]
[injection_rate_massfrac_H2O]
initial_condition = 0.994175112
[]
[rate_H2O]
[]
[rate_Na]
[]
[rate_Cl]
[]
[rate_SiO2]
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosityConst # this simulation has no porosity changes from dissolution
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '0 0 0 0 0 0 0 0 0'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[]
[Preconditioning]
active = typically_efficient
[typically_efficient]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = ' hypre boomeramg'
[]
[strong]
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 ilu NONZERO 2'
[]
[probably_too_strong]
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 = 7.76E6 # 90 days
dt = 1E5
[]
[Outputs]
exodus = true
[]
[MultiApps]
[react]
type = TransientMultiApp
input_files = aquifer_geochemistry.i
clone_master_mesh = true
execute_on = 'timestep_end'
[]
[]
[Transfers]
[changes_due_to_flow]
type = MultiAppCopyTransfer
direction = to_multiapp
source_variable = 'rate_H2O rate_Na rate_Cl rate_SiO2 temperature'
variable = 'pf_rate_H2O pf_rate_Na pf_rate_Cl pf_rate_SiO2 temperature'
multi_app = react
[]
[massfrac_from_geochem]
type = MultiAppCopyTransfer
direction = from_multiapp
source_variable = 'massfrac_Na massfrac_Cl massfrac_SiO2'
variable = 'f0 f1 f2'
multi_app = react
[]
[]
(modules/combined/examples/geochem-porous_flow/geotes_weber_tensleep/porous_flow.i)
#########################################
# #
# File written by create_input_files.py #
# #
#########################################
# PorousFlow simulation of injection and production in a simplified GeoTES aquifer
# Much of this file is standard porous-flow stuff. The unusual aspects are:
# - transfer of the rates of changes of each species (kg.s) to the aquifer_geochemistry.i simulation. This is achieved by saving these changes from the PorousFlowMassTimeDerivative residuals
# - transfer of the temperature field to the aquifer_geochemistry.i simulation
# Interesting behaviour can be simulated by this file without its 'parent' simulation, exchanger.i. exchanger.i provides mass-fractions injected via the injection_rate_massfrac_* variables, but since these are more-or-less constant throughout the duration of the exchanger.i simulation, the initial_conditions specified below may be used. Similar, exchanger.i provides injection_temperature, but that is also constant.
injection_rate = -0.02 # kg/s/m, negative because injection as a source
production_rate = 0.02 # kg/s/m, this is about the maximum that can be sustained by the aquifer, with its fairly low permeability, without porepressure becoming negative
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 3
xmin = -75
xmax = 75
ymin = 0
ymax = 40
zmin = -25
zmax = 25
nx = 15
ny = 4
nz = 5
[]
[aquifer]
type = ParsedSubdomainMeshGenerator
input = gen
block_id = 1
block_name = aquifer
combinatorial_geometry = 'z >= -5 & z <= 5'
[]
[injection_nodes]
input = aquifer
type = ExtraNodesetGenerator
new_boundary = injection_nodes
coord = '-25 0 -5; -25 0 5'
[]
[production_nodes]
input = injection_nodes
type = ExtraNodesetGenerator
new_boundary = production_nodes
coord = '25 0 -5; 25 0 5'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 -10'
[]
[BCs]
[injection_temperature]
type = MatchedValueBC
variable = temperature
v = injection_temperature
boundary = injection_nodes
[]
[]
[Modules]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0
bulk_modulus = 2E9
viscosity = 1E-3
density0 = 1000
cv = 4000.0
cp = 4000.0
[]
[]
[]
[PorousFlowFullySaturated]
coupling_type = ThermoHydro
porepressure = porepressure
temperature = temperature
mass_fraction_vars = 'f_H f_Cl f_SO4 f_HCO3 f_SiO2aq f_Al f_Ca f_Mg f_Fe f_K f_Na f_Sr f_F f_BOH f_Br f_Ba f_Li f_NO3 f_O2aq '
save_component_rate_in = 'rate_H rate_Cl rate_SO4 rate_HCO3 rate_SiO2aq rate_Al rate_Ca rate_Mg rate_Fe rate_K rate_Na rate_Sr rate_F rate_BOH rate_Br rate_Ba rate_Li rate_NO3 rate_O2aq rate_H2O' # change in kg at every node / dt
fp = the_simple_fluid
temperature_unit = Celsius
[]
[Materials]
[porosity_caps]
type = PorousFlowPorosityConst # this simulation has no porosity changes from dissolution
block = 0
porosity = 0.01
[]
[porosity_aquifer]
type = PorousFlowPorosityConst # this simulation has no porosity changes from dissolution
block = aquifer
porosity = 0.063
[]
[permeability_caps]
type = PorousFlowPermeabilityConst
block = 0
permeability = '1E-18 0 0 0 1E-18 0 0 0 1E-18'
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1.7E-15 0 0 0 1.7E-15 0 0 0 4.1E-16'
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '0 0 0 0 0 0 0 0 0'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[]
[Preconditioning]
active = typically_efficient
[typically_efficient]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = ' hypre boomeramg'
[]
[strong]
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 ilu NONZERO 2'
[]
[probably_too_strong]
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 = 7.76E6 # 90 days
[TimeStepper]
type = FunctionDT
function = 'min(3E4, max(1E4, 0.2 * t))'
[]
[]
[Outputs]
exodus = true
[]
[Variables]
[f_H]
initial_condition = -2.952985071156e-06
[]
[f_Cl]
initial_condition = 0.04870664551708
[]
[f_SO4]
initial_condition = 0.0060359986852517
[]
[f_HCO3]
initial_condition = 5.0897287594019e-05
[]
[f_SiO2aq]
initial_condition = 3.0246609868421e-05
[]
[f_Al]
initial_condition = 3.268028901929e-08
[]
[f_Ca]
initial_condition = 0.00082159428184586
[]
[f_Mg]
initial_condition = 1.8546347062146e-05
[]
[f_Fe]
initial_condition = 4.3291908204093e-05
[]
[f_K]
initial_condition = 6.8434768308898e-05
[]
[f_Na]
initial_condition = 0.033298053919671
[]
[f_Sr]
initial_condition = 1.2771866652177e-05
[]
[f_F]
initial_condition = 5.5648860174073e-06
[]
[f_BOH]
initial_condition = 0.0003758574621917
[]
[f_Br]
initial_condition = 9.0315286107068e-05
[]
[f_Ba]
initial_condition = 1.5637460875161e-07
[]
[f_Li]
initial_condition = 8.3017067912701e-05
[]
[f_NO3]
initial_condition = 0.00010958455036169
[]
[f_O2aq]
initial_condition = -7.0806852373351e-05
[]
[porepressure]
initial_condition = 30E6
[]
[temperature]
initial_condition = 92
scaling = 1E-6 # fluid enthalpy is roughly 1E6
[]
[]
[DiracKernels]
[inject_H]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_H
point_file = injection.bh
variable = f_H
[]
[inject_Cl]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Cl
point_file = injection.bh
variable = f_Cl
[]
[inject_SO4]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_SO4
point_file = injection.bh
variable = f_SO4
[]
[inject_HCO3]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_HCO3
point_file = injection.bh
variable = f_HCO3
[]
[inject_SiO2aq]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_SiO2aq
point_file = injection.bh
variable = f_SiO2aq
[]
[inject_Al]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Al
point_file = injection.bh
variable = f_Al
[]
[inject_Ca]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Ca
point_file = injection.bh
variable = f_Ca
[]
[inject_Mg]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Mg
point_file = injection.bh
variable = f_Mg
[]
[inject_Fe]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Fe
point_file = injection.bh
variable = f_Fe
[]
[inject_K]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_K
point_file = injection.bh
variable = f_K
[]
[inject_Na]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Na
point_file = injection.bh
variable = f_Na
[]
[inject_Sr]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Sr
point_file = injection.bh
variable = f_Sr
[]
[inject_F]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_F
point_file = injection.bh
variable = f_F
[]
[inject_BOH]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_BOH
point_file = injection.bh
variable = f_BOH
[]
[inject_Br]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Br
point_file = injection.bh
variable = f_Br
[]
[inject_Ba]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Ba
point_file = injection.bh
variable = f_Ba
[]
[inject_Li]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_Li
point_file = injection.bh
variable = f_Li
[]
[inject_NO3]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_NO3
point_file = injection.bh
variable = f_NO3
[]
[inject_O2aq]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_O2aq
point_file = injection.bh
variable = f_O2aq
[]
[inject_H2O]
type = PorousFlowPolyLineSink
SumQuantityUO = injected_mass
fluxes = ${injection_rate}
p_or_t_vals = 0.0
multiplying_var = injection_rate_massfrac_H2O
point_file = injection.bh
variable = porepressure
[]
[produce_H]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_H
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 0
point_file = production.bh
variable = f_H
[]
[produce_Cl]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Cl
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 1
point_file = production.bh
variable = f_Cl
[]
[produce_SO4]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_SO4
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 2
point_file = production.bh
variable = f_SO4
[]
[produce_HCO3]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_HCO3
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 3
point_file = production.bh
variable = f_HCO3
[]
[produce_SiO2aq]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_SiO2aq
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 4
point_file = production.bh
variable = f_SiO2aq
[]
[produce_Al]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Al
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 5
point_file = production.bh
variable = f_Al
[]
[produce_Ca]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Ca
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 6
point_file = production.bh
variable = f_Ca
[]
[produce_Mg]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Mg
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 7
point_file = production.bh
variable = f_Mg
[]
[produce_Fe]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Fe
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 8
point_file = production.bh
variable = f_Fe
[]
[produce_K]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_K
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 9
point_file = production.bh
variable = f_K
[]
[produce_Na]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Na
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 10
point_file = production.bh
variable = f_Na
[]
[produce_Sr]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Sr
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 11
point_file = production.bh
variable = f_Sr
[]
[produce_F]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_F
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 12
point_file = production.bh
variable = f_F
[]
[produce_BOH]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_BOH
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 13
point_file = production.bh
variable = f_BOH
[]
[produce_Br]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Br
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 14
point_file = production.bh
variable = f_Br
[]
[produce_Ba]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Ba
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 15
point_file = production.bh
variable = f_Ba
[]
[produce_Li]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_Li
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 16
point_file = production.bh
variable = f_Li
[]
[produce_NO3]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_NO3
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 17
point_file = production.bh
variable = f_NO3
[]
[produce_O2aq]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_O2aq
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 18
point_file = production.bh
variable = f_O2aq
[]
[produce_H2O]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_mass_H2O
fluxes = ${production_rate}
p_or_t_vals = 0.0
mass_fraction_component = 19
point_file = production.bh
variable = porepressure
[]
[produce_heat]
type = PorousFlowPolyLineSink
SumQuantityUO = produced_heat
fluxes = ${production_rate}
p_or_t_vals = 0.0
use_enthalpy = true
point_file = production.bh
variable = temperature
[]
[]
[UserObjects]
[injected_mass]
type = PorousFlowSumQuantity
[]
[produced_mass_H]
type = PorousFlowSumQuantity
[]
[produced_mass_Cl]
type = PorousFlowSumQuantity
[]
[produced_mass_SO4]
type = PorousFlowSumQuantity
[]
[produced_mass_HCO3]
type = PorousFlowSumQuantity
[]
[produced_mass_SiO2aq]
type = PorousFlowSumQuantity
[]
[produced_mass_Al]
type = PorousFlowSumQuantity
[]
[produced_mass_Ca]
type = PorousFlowSumQuantity
[]
[produced_mass_Mg]
type = PorousFlowSumQuantity
[]
[produced_mass_Fe]
type = PorousFlowSumQuantity
[]
[produced_mass_K]
type = PorousFlowSumQuantity
[]
[produced_mass_Na]
type = PorousFlowSumQuantity
[]
[produced_mass_Sr]
type = PorousFlowSumQuantity
[]
[produced_mass_F]
type = PorousFlowSumQuantity
[]
[produced_mass_BOH]
type = PorousFlowSumQuantity
[]
[produced_mass_Br]
type = PorousFlowSumQuantity
[]
[produced_mass_Ba]
type = PorousFlowSumQuantity
[]
[produced_mass_Li]
type = PorousFlowSumQuantity
[]
[produced_mass_NO3]
type = PorousFlowSumQuantity
[]
[produced_mass_O2aq]
type = PorousFlowSumQuantity
[]
[produced_mass_H2O]
type = PorousFlowSumQuantity
[]
[produced_heat]
type = PorousFlowSumQuantity
[]
[]
[Postprocessors]
[dt]
type = TimestepSize
execute_on = TIMESTEP_BEGIN
[]
[tot_kg_injected_this_timestep]
type = PorousFlowPlotQuantity
uo = injected_mass
[]
[kg_H_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_H
[]
[kg_Cl_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Cl
[]
[kg_SO4_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_SO4
[]
[kg_HCO3_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_HCO3
[]
[kg_SiO2aq_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_SiO2aq
[]
[kg_Al_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Al
[]
[kg_Ca_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Ca
[]
[kg_Mg_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Mg
[]
[kg_Fe_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Fe
[]
[kg_K_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_K
[]
[kg_Na_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Na
[]
[kg_Sr_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Sr
[]
[kg_F_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_F
[]
[kg_BOH_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_BOH
[]
[kg_Br_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Br
[]
[kg_Ba_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Ba
[]
[kg_Li_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_Li
[]
[kg_NO3_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_NO3
[]
[kg_O2aq_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_O2aq
[]
[kg_H2O_produced_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_mass_H2O
[]
[mole_rate_H_produced]
type = FunctionValuePostprocessor
function = moles_H
[]
[mole_rate_Cl_produced]
type = FunctionValuePostprocessor
function = moles_Cl
[]
[mole_rate_SO4_produced]
type = FunctionValuePostprocessor
function = moles_SO4
[]
[mole_rate_HCO3_produced]
type = FunctionValuePostprocessor
function = moles_HCO3
[]
[mole_rate_SiO2aq_produced]
type = FunctionValuePostprocessor
function = moles_SiO2aq
[]
[mole_rate_Al_produced]
type = FunctionValuePostprocessor
function = moles_Al
[]
[mole_rate_Ca_produced]
type = FunctionValuePostprocessor
function = moles_Ca
[]
[mole_rate_Mg_produced]
type = FunctionValuePostprocessor
function = moles_Mg
[]
[mole_rate_Fe_produced]
type = FunctionValuePostprocessor
function = moles_Fe
[]
[mole_rate_K_produced]
type = FunctionValuePostprocessor
function = moles_K
[]
[mole_rate_Na_produced]
type = FunctionValuePostprocessor
function = moles_Na
[]
[mole_rate_Sr_produced]
type = FunctionValuePostprocessor
function = moles_Sr
[]
[mole_rate_F_produced]
type = FunctionValuePostprocessor
function = moles_F
[]
[mole_rate_BOH_produced]
type = FunctionValuePostprocessor
function = moles_BOH
[]
[mole_rate_Br_produced]
type = FunctionValuePostprocessor
function = moles_Br
[]
[mole_rate_Ba_produced]
type = FunctionValuePostprocessor
function = moles_Ba
[]
[mole_rate_Li_produced]
type = FunctionValuePostprocessor
function = moles_Li
[]
[mole_rate_NO3_produced]
type = FunctionValuePostprocessor
function = moles_NO3
[]
[mole_rate_O2aq_produced]
type = FunctionValuePostprocessor
function = moles_O2aq
[]
[mole_rate_H2O_produced]
type = FunctionValuePostprocessor
function = moles_H2O
[]
[heat_joules_extracted_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_heat
[]
[production_temperature]
type = AverageNodalVariableValue
boundary = production_nodes
variable = temperature
[]
[]
[Functions]
[moles_H]
type = ParsedFunction
vars = 'kg_H dt'
vals = 'kg_H_produced_this_timestep dt'
value = 'kg_H * 1000 / 1.0079 / dt'
[]
[moles_Cl]
type = ParsedFunction
vars = 'kg_Cl dt'
vals = 'kg_Cl_produced_this_timestep dt'
value = 'kg_Cl * 1000 / 35.453 / dt'
[]
[moles_SO4]
type = ParsedFunction
vars = 'kg_SO4 dt'
vals = 'kg_SO4_produced_this_timestep dt'
value = 'kg_SO4 * 1000 / 96.0576 / dt'
[]
[moles_HCO3]
type = ParsedFunction
vars = 'kg_HCO3 dt'
vals = 'kg_HCO3_produced_this_timestep dt'
value = 'kg_HCO3 * 1000 / 61.0171 / dt'
[]
[moles_SiO2aq]
type = ParsedFunction
vars = 'kg_SiO2aq dt'
vals = 'kg_SiO2aq_produced_this_timestep dt'
value = 'kg_SiO2aq * 1000 / 60.0843 / dt'
[]
[moles_Al]
type = ParsedFunction
vars = 'kg_Al dt'
vals = 'kg_Al_produced_this_timestep dt'
value = 'kg_Al * 1000 / 26.9815 / dt'
[]
[moles_Ca]
type = ParsedFunction
vars = 'kg_Ca dt'
vals = 'kg_Ca_produced_this_timestep dt'
value = 'kg_Ca * 1000 / 40.08 / dt'
[]
[moles_Mg]
type = ParsedFunction
vars = 'kg_Mg dt'
vals = 'kg_Mg_produced_this_timestep dt'
value = 'kg_Mg * 1000 / 24.305 / dt'
[]
[moles_Fe]
type = ParsedFunction
vars = 'kg_Fe dt'
vals = 'kg_Fe_produced_this_timestep dt'
value = 'kg_Fe * 1000 / 55.847 / dt'
[]
[moles_K]
type = ParsedFunction
vars = 'kg_K dt'
vals = 'kg_K_produced_this_timestep dt'
value = 'kg_K * 1000 / 39.0983 / dt'
[]
[moles_Na]
type = ParsedFunction
vars = 'kg_Na dt'
vals = 'kg_Na_produced_this_timestep dt'
value = 'kg_Na * 1000 / 22.9898 / dt'
[]
[moles_Sr]
type = ParsedFunction
vars = 'kg_Sr dt'
vals = 'kg_Sr_produced_this_timestep dt'
value = 'kg_Sr * 1000 / 87.62 / dt'
[]
[moles_F]
type = ParsedFunction
vars = 'kg_F dt'
vals = 'kg_F_produced_this_timestep dt'
value = 'kg_F * 1000 / 18.9984 / dt'
[]
[moles_BOH]
type = ParsedFunction
vars = 'kg_BOH dt'
vals = 'kg_BOH_produced_this_timestep dt'
value = 'kg_BOH * 1000 / 61.8329 / dt'
[]
[moles_Br]
type = ParsedFunction
vars = 'kg_Br dt'
vals = 'kg_Br_produced_this_timestep dt'
value = 'kg_Br * 1000 / 79.904 / dt'
[]
[moles_Ba]
type = ParsedFunction
vars = 'kg_Ba dt'
vals = 'kg_Ba_produced_this_timestep dt'
value = 'kg_Ba * 1000 / 137.33 / dt'
[]
[moles_Li]
type = ParsedFunction
vars = 'kg_Li dt'
vals = 'kg_Li_produced_this_timestep dt'
value = 'kg_Li * 1000 / 6.941 / dt'
[]
[moles_NO3]
type = ParsedFunction
vars = 'kg_NO3 dt'
vals = 'kg_NO3_produced_this_timestep dt'
value = 'kg_NO3 * 1000 / 62.0049 / dt'
[]
[moles_O2aq]
type = ParsedFunction
vars = 'kg_O2aq dt'
vals = 'kg_O2aq_produced_this_timestep dt'
value = 'kg_O2aq * 1000 / 31.9988 / dt'
[]
[moles_H2O]
type = ParsedFunction
vars = 'kg_H2O dt'
vals = 'kg_H2O_produced_this_timestep dt'
value = 'kg_H2O * 1000 / 18.01801802 / dt'
[]
[]
[AuxVariables]
[injection_temperature]
initial_condition = 92
[]
[injection_rate_massfrac_H]
initial_condition = -2.952985071156e-06
[]
[injection_rate_massfrac_Cl]
initial_condition = 0.04870664551708
[]
[injection_rate_massfrac_SO4]
initial_condition = 0.0060359986852517
[]
[injection_rate_massfrac_HCO3]
initial_condition = 5.0897287594019e-05
[]
[injection_rate_massfrac_SiO2aq]
initial_condition = 3.0246609868421e-05
[]
[injection_rate_massfrac_Al]
initial_condition = 3.268028901929e-08
[]
[injection_rate_massfrac_Ca]
initial_condition = 0.00082159428184586
[]
[injection_rate_massfrac_Mg]
initial_condition = 1.8546347062146e-05
[]
[injection_rate_massfrac_Fe]
initial_condition = 4.3291908204093e-05
[]
[injection_rate_massfrac_K]
initial_condition = 6.8434768308898e-05
[]
[injection_rate_massfrac_Na]
initial_condition = 0.033298053919671
[]
[injection_rate_massfrac_Sr]
initial_condition = 1.2771866652177e-05
[]
[injection_rate_massfrac_F]
initial_condition = 5.5648860174073e-06
[]
[injection_rate_massfrac_BOH]
initial_condition = 0.0003758574621917
[]
[injection_rate_massfrac_Br]
initial_condition = 9.0315286107068e-05
[]
[injection_rate_massfrac_Ba]
initial_condition = 1.5637460875161e-07
[]
[injection_rate_massfrac_Li]
initial_condition = 8.3017067912701e-05
[]
[injection_rate_massfrac_NO3]
initial_condition = 0.00010958455036169
[]
[injection_rate_massfrac_O2aq]
initial_condition = -7.0806852373351e-05
[]
[injection_rate_massfrac_H2O]
initial_condition = 0.91032275033842
[]
[rate_H]
[]
[rate_Cl]
[]
[rate_SO4]
[]
[rate_HCO3]
[]
[rate_SiO2aq]
[]
[rate_Al]
[]
[rate_Ca]
[]
[rate_Mg]
[]
[rate_Fe]
[]
[rate_K]
[]
[rate_Na]
[]
[rate_Sr]
[]
[rate_F]
[]
[rate_BOH]
[]
[rate_Br]
[]
[rate_Ba]
[]
[rate_Li]
[]
[rate_NO3]
[]
[rate_O2aq]
[]
[rate_H2O]
[]
[]
[MultiApps]
[react]
type = TransientMultiApp
input_files = aquifer_geochemistry.i
clone_master_mesh = true
execute_on = 'timestep_end'
[]
[]
[Transfers]
[changes_due_to_flow]
type = MultiAppCopyTransfer
direction = to_multiapp
source_variable = 'rate_H rate_Cl rate_SO4 rate_HCO3 rate_SiO2aq rate_Al rate_Ca rate_Mg rate_Fe rate_K rate_Na rate_Sr rate_F rate_BOH rate_Br rate_Ba rate_Li rate_NO3 rate_O2aq rate_H2O temperature'
variable = 'pf_rate_H pf_rate_Cl pf_rate_SO4 pf_rate_HCO3 pf_rate_SiO2aq pf_rate_Al pf_rate_Ca pf_rate_Mg pf_rate_Fe pf_rate_K pf_rate_Na pf_rate_Sr pf_rate_F pf_rate_BOH pf_rate_Br pf_rate_Ba pf_rate_Li pf_rate_NO3 pf_rate_O2aq pf_rate_H2O temperature'
multi_app = react
[]
[massfrac_from_geochem]
type = MultiAppCopyTransfer
direction = from_multiapp
source_variable = 'massfrac_H massfrac_Cl massfrac_SO4 massfrac_HCO3 massfrac_SiO2aq massfrac_Al massfrac_Ca massfrac_Mg massfrac_Fe massfrac_K massfrac_Na massfrac_Sr massfrac_F massfrac_BOH massfrac_Br massfrac_Ba massfrac_Li massfrac_NO3 massfrac_O2aq '
variable = 'f_H f_Cl f_SO4 f_HCO3 f_SiO2aq f_Al f_Ca f_Mg f_Fe f_K f_Na f_Sr f_F f_BOH f_Br f_Ba f_Li f_NO3 f_O2aq '
multi_app = react
[]
[]
(modules/porous_flow/test/tests/dirackernels/bh_except15.i)
# fully-saturated
# production
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/examples/groundwater/ex01.i)
# Groundwater extraction example.
# System consists of two confined aquifers separated by an aquitard
# There is a hydraulic gradient in the upper aquifer
# A well extracts water from the lower aquifer, and the impact on the upper aquifer is observed
# In the center of the model, the roof of the upper aquifer sits 70m below the local water table
[Mesh]
[basic_mesh]
type = GeneratedMeshGenerator
dim = 3
xmin = -50
xmax = 50
nx = 20
ymin = -25
ymax = 25
ny = 10
zmin = -100
zmax = -70
nz = 3
[]
[lower_aquifer]
type = SubdomainBoundingBoxGenerator
input = basic_mesh
block_id = 1
block_name = lower_aquifer
bottom_left = '-1000 -500 -100'
top_right = '1000 500 -90'
[]
[aquitard]
type = SubdomainBoundingBoxGenerator
input = lower_aquifer
block_id = 2
block_name = aquitard
bottom_left = '-1000 -500 -90'
top_right = '1000 500 -80'
[]
[upper_aquifer]
type = SubdomainBoundingBoxGenerator
input = aquitard
block_id = 3
block_name = upper_aquifer
bottom_left = '-1000 -500 -80'
top_right = '1000 500 -70'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
function = insitu_pp
[]
[]
[BCs]
[pp]
type = FunctionDirichletBC
variable = pp
function = insitu_pp
boundary = 'left right top bottom front back'
[]
[]
[Functions]
[upper_aquifer_head]
type = ParsedFunction
value = '10 + x / 200'
[]
[lower_aquifer_head]
type = ParsedFunction
value = '20'
[]
[insitu_head]
type = ParsedFunction
vals = 'lower_aquifer_head upper_aquifer_head'
vars = 'low up'
value = 'if(z <= -90, low, if(z >= -80, up, (up * (z + 90) - low * (z + 80)) / (10.0)))'
[]
[insitu_pp]
type = ParsedFunction
vals = 'insitu_head'
vars = 'h'
value = '(h - z) * 1E4'
[]
[l_rate]
type = ParsedFunction
vals = 'm3_produced dt'
vars = 'm3_produced dt'
value = '1000 * m3_produced / dt'
[]
[]
[AuxVariables]
[insitu_head]
[]
[head_change]
[]
[]
[AuxKernels]
[insitu_head]
type = FunctionAux
variable = insitu_head
function = insitu_head
[]
[head_change]
type = ParsedAux
args = 'pp insitu_head'
use_xyzt = true
function = 'pp / 1E4 + z - insitu_head'
variable = head_change
[]
[]
[Postprocessors]
[m3_produced]
type = PorousFlowPlotQuantity
uo = volume_extracted
outputs = 'none'
[]
[dt]
type = TimestepSize
outputs = 'none'
[]
[l_per_s]
type = FunctionValuePostprocessor
function = l_rate
[]
[]
[VectorPostprocessors]
[drawdown]
type = LineValueSampler
variable = head_change
start_point = '-50 0 -75'
end_point = '50 0 -75'
num_points = 101
sort_by = x
[]
[]
[PorousFlowBasicTHM]
fp = simple_fluid
gravity = '0 0 -10'
porepressure = pp
multiply_by_density = false
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
# the following mean that density = 1000 * exp(P / 1E15) ~ 1000
thermal_expansion = 0
bulk_modulus = 1E15
[]
[]
[]
[Materials]
[porosity_aquifers]
type = PorousFlowPorosityConst
porosity = 0.05
block = 'upper_aquifer lower_aquifer'
[]
[porosity_aquitard]
type = PorousFlowPorosityConst
porosity = 0.2
block = aquitard
[]
[biot_mod]
type = PorousFlowConstantBiotModulus
fluid_bulk_modulus = 2E9
biot_coefficient = 1.0
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
block = 'upper_aquifer lower_aquifer'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst
permeability = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
block = aquitard
[]
[]
[DiracKernels]
[sink]
type = PorousFlowPolyLineSink
SumQuantityUO = volume_extracted
point_file = ex01.bh_lower
line_length = 10
variable = pp
# following produces a flux of 0 m^3(water)/m(borehole length)/s if porepressure = 0, and a flux of 1 m^3/m/s if porepressure = 1E9
p_or_t_vals = '0 1E9'
fluxes = '0 1'
[]
[]
[UserObjects]
[volume_extracted]
type = PorousFlowSumQuantity
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
[TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 1.1E5
[]
end_time = 3.456E5 # 4 days
nl_abs_tol = 1E-13
[]
[Outputs]
[csv]
type = CSV
file_base = ex01_lower_extraction
execute_on = final
[]
[]
(modules/porous_flow/test/tests/dirackernels/pls03_action.i)
# Test that the upwinding works correctly.
#
# A poly-line sink sits at the centre of the element.
# It has length=4 and weight=0.5, and extracts fluid
# at a constant rate of
# (1 * relative_permeability) kg.m^-1.s^-1
# Since it sits at the centre of the element, it extracts
# equally from each node, so the rate of extraction from
# each node is
# (0.5 * relative_permeability) kg.s^-1
# including the length and weight effects.
#
# There is no fluid flow.
#
# The initial conditions are such that all nodes have
# relative_permeability=0, except for one which has
# relative_permeaility = 1. Therefore, all nodes should
# remain at their initial porepressure, except the one.
#
# The porosity is 0.1, and the elemental volume is 2,
# so the fluid mass at the node in question = 0.2 * density / 4,
# where the 4 is the number of nodes in the element.
# In this simulation density = dens0 * exp(P / bulk), with
# dens0 = 100, and bulk = 20 MPa.
# The initial porepressure P0 = 10 MPa, so the final (after
# 1 second of simulation) is
# P(t=1) = 8.748592 MPa
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = 0
xmax = 2
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[Modules]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 2.0E7
viscosity = 1.0
density0 = 100.0
[]
[]
[]
[PorousFlowUnsaturated]
porepressure = pp
gravity = '0 0 0'
fp = the_simple_fluid
van_genuchten_alpha = 1.0E-7
van_genuchten_m = 0.5
relative_permeability_exponent = 2
residual_saturation = 0.99
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
#function = if((x<1)&(y<0.5),1E7,-1E7)
function = if((x<1)&(y>0.5),1E7,-1E7)
#function = if((x>1)&(y<0.5),1E7,-1E7)
#function = if((x>1)&(y>0.5),1E7,-1E7)
[]
[]
[UserObjects]
[pls_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[]
[Materials]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '0 0 0 0 0 0 0 0 0'
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[DiracKernels]
[pls]
type = PorousFlowPolyLineSink
fluid_phase = 0
point_file = pls03.bh
use_relative_permeability = true
line_length = 4
SumQuantityUO = pls_total_outflow_mass
variable = pp
p_or_t_vals = '0 1E7'
fluxes = '1 1'
[]
[]
[Postprocessors]
[pls_report]
type = PorousFlowPlotQuantity
uo = pls_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p00]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[p01]
type = PointValue
variable = pp
point = '0 1 0'
execute_on = timestep_end
[]
[p20]
type = PointValue
variable = pp
point = '2 0 0'
execute_on = timestep_end
[]
[p21]
type = PointValue
variable = pp
point = '2 1 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 pls_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 1
dt = 1
solve_type = NEWTON
[]
[Outputs]
file_base = pls03_action
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/bh_except03.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
at_nodes = true # Needed to force expected error
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except07.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
use_mobility = true
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except13.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = coincident_points.bh
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/examples/groundwater/ex02_steady_state.i)
# Steady-state groundwater model. See groundwater_models.md for a detailed description
[Mesh]
[basic_mesh]
# mesh create by external program: lies within -500<=x<=500 and -200<=y<=200, with varying z
type = FileMeshGenerator
file = ex02_mesh.e
[]
[name_blocks]
type = RenameBlockGenerator
input = basic_mesh
old_block_id = '2 3 4'
new_block_name = 'bot_aquifer aquitard top_aquifer'
[]
[zmax]
type = SideSetsFromNormalsGenerator
input = name_blocks
new_boundary = zmax
normals = '0 0 1'
[]
[xmin_bot_aquifer]
type = ParsedGenerateSideset
input = zmax
included_subdomain_ids = 2
normal = '-1 0 0'
combinatorial_geometry = 'x <= -500.0'
new_sideset_name = xmin_bot_aquifer
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
function = initial_pp
[]
[]
[BCs]
[rainfall_recharge]
type = PorousFlowSink
boundary = zmax
variable = pp
flux_function = -1E-6 # recharge of 0.1mm/day = 1E-4m3/m2/day = 0.1kg/m2/day ~ 1E-6kg/m2/s
[]
[evapotranspiration]
type = PorousFlowHalfCubicSink
boundary = zmax
variable = pp
center = 0.0
cutoff = -5E4 # roots of depth 5m. 5m of water = 5E4 Pa
use_mobility = true
fluid_phase = 0
# Assume pan evaporation of 4mm/day = 4E-3m3/m2/day = 4kg/m2/day ~ 4E-5kg/m2/s
# Assume that if permeability was 1E-10m^2 and water table at topography then ET acts as pan strength
# Because use_mobility = true, then 4E-5 = maximum_flux = max * perm * density / visc = max * 1E-4, so max = 40
max = 40
[]
[]
[DiracKernels]
[river]
type = PorousFlowPolyLineSink
SumQuantityUO = baseflow
point_file = ex02_river.bh
# Assume a perennial river.
# Assume the river has an incision depth of 1m and a stage height of 1.5m, and these are constant in time and uniform over the whole model. Hence, if groundwater head is 0.5m (5000Pa) there will be no baseflow and leakage.
p_or_t_vals = '-999995000 5000 1000005000'
# Assume the riverbed conductance, k_zz*density*river_segment_length*river_width/riverbed_thickness/viscosity = 1E-6*river_segment_length kg/Pa/s
fluxes = '-1E3 0 1E3'
variable = pp
[]
[]
[Functions]
[initial_pp]
type = SolutionFunction
scale_factor = 1E4
from_variable = cosflow_depth
solution = initial_mesh
[]
[baseflow_rate]
type = ParsedFunction
vars = 'baseflow_kg dt'
vals = 'baseflow_kg dt'
value = 'baseflow_kg / dt * 24.0 * 3600.0 / 400.0'
[]
[]
[PorousFlowUnsaturated]
fp = simple_fluid
porepressure = pp
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[]
[Materials]
[porosity_everywhere]
type = PorousFlowPorosityConst
porosity = 0.05
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst
block = 'top_aquifer bot_aquifer'
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-13'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst
block = aquitard
permeability = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
[]
[]
[UserObjects]
[initial_mesh]
type = SolutionUserObject
execute_on = INITIAL
mesh = ex02_mesh.e
timestep = LATEST
system_variables = cosflow_depth
[]
[baseflow]
type = PorousFlowSumQuantity
[]
[]
[Postprocessors]
[baseflow_kg]
type = PorousFlowPlotQuantity
uo = baseflow
outputs = 'none'
[]
[dt]
type = TimestepSize
outputs = 'none'
[]
[baseflow_l_per_m_per_day]
type = FunctionValuePostprocessor
function = baseflow_rate
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
# following 2 lines are not mandatory, but illustrate a popular preconditioner choice in groundwater models
petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap'
petsc_options_value = ' asm ilu 2 '
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1E6
[TimeStepper]
type = FunctionDT
function = 'max(1E6, t)'
[]
end_time = 1E12
nl_abs_tol = 1E-13
[]
[Outputs]
print_linear_residuals = false
[ex]
type = Exodus
execute_on = final
[]
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/dirackernels/bh04.i)
# fully-saturated
# production
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Functions]
[dts]
type = PiecewiseLinear
y = '1E-2 1E-1 1 1E1 1E2 1E3'
x = '0 1E-1 1 1E1 1E2 1E3'
[]
[]
[Variables]
[pp]
initial_condition = 0
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityFLAC
m = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
variable = pp
SumQuantityUO = borehole_total_outflow_mass
point_file = bh02.bh
fluid_phase = 0
bottom_p_or_t = -1E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 1E3
solve_type = NEWTON
[TimeStepper]
type = FunctionDT
function = dts
[]
[]
[Outputs]
file_base = bh04
exodus = false
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/dirackernels/injection_production.i)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 3
nx = 10
ny = 10
nz = 1
xmin = -50
xmax = 50
ymin = -50
ymax = 50
zmin = 0
zmax = 10
[]
[central_nodes]
input = gen
type = ExtraNodesetGenerator
new_boundary = central_nodes
coord = '0 0 0; 0 0 10'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 20E6
[]
[temperature]
initial_condition = 400
scaling = 1E-6 # fluid enthalpy is roughly 1E6
[]
[]
[BCs]
[injection_temperature]
type = DirichletBC
variable = temperature
value = 300
boundary = central_nodes
[]
[]
[DiracKernels]
[fluid_injection]
type = PorousFlowPeacemanBorehole
variable = porepressure
SumQuantityUO = injected_mass
point_file = injection.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 21E6
unit_weight = '0 0 0'
use_mobility = true
character = -1
[]
[fluid_production]
type = PorousFlowPeacemanBorehole
variable = porepressure
SumQuantityUO = produced_mass
point_file = production.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[remove_heat_at_production_well]
type = PorousFlowPeacemanBorehole
variable = temperature
SumQuantityUO = produced_heat
point_file = production.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
use_enthalpy = true
character = 1
[]
[]
[UserObjects]
[injected_mass]
type = PorousFlowSumQuantity
[]
[produced_mass]
type = PorousFlowSumQuantity
[]
[produced_heat]
type = PorousFlowSumQuantity
[]
[]
[Postprocessors]
[heat_joules_extracted_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_heat
[]
[]
[Modules]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 2E-4
bulk_modulus = 2E9
viscosity = 1E-3
density0 = 1000
cv = 4000.0
cp = 4000.0
[]
[]
[]
[PorousFlowUnsaturated]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 1E-10
fluid_bulk_modulus = 2E9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
fluid_coefficient = 5E-6
drained_coefficient = 2E-4
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1 0 0 0 1 0 0 0 1'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 2E6
dt = 2E5
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/dirackernels/bh_except09.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
compute_enthalpy = false
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
use_mobility = true
use_enthalpy = true
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]
(modules/porous_flow/test/tests/dirackernels/bh_except06.i)
# PorousFlowPeacemanBorehole exception test
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 1E7
[]
[]
[Kernels]
[mass0]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[Modules]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
bottom_p_or_t = 0
fluid_phase = 0
point_file = bh02.bh
use_mobility = true
SumQuantityUO = borehole_total_outflow_mass
variable = pp
unit_weight = '0 0 0'
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[fluid_mass0]
type = PorousFlowFluidMass
execute_on = timestep_begin
[]
[fluid_mass1]
type = PorousFlowFluidMass
execute_on = timestep_end
[]
[zmass_error]
type = FunctionValuePostprocessor
function = mass_bal_fcn
execute_on = timestep_end
[]
[p0]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = timestep_end
[]
[]
[Functions]
[mass_bal_fcn]
type = ParsedFunction
value = abs((a-c+d)/2/(a+c))
vars = 'a c d'
vals = 'fluid_mass1 fluid_mass0 bh_report'
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1E-10 1E-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 1E-2
solve_type = NEWTON
[]