- PorousFlowDictatorThe UserObject that holds the list of PorousFlow variable names
C++ Type:UserObjectName
Controllable:No
Description:The UserObject that holds the list of PorousFlow variable names
- SumQuantityUOUser Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step.
C++ Type:UserObjectName
Controllable:No
Description:User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step.
- T_inThe incoming fluid temperature
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:The incoming fluid temperature
- fluxesTuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed). These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1. A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature). If a quad-point pressure is less than the first pressure value, the first flux value is used. If quad-point pressure exceeds the final pressure value, the final flux value is used. This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE.
C++ Type:std::vector<Real>
Unit:(no unit assumed)
Controllable:No
Description:Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed). These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1. A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature). If a quad-point pressure is less than the first pressure value, the first flux value is used. If quad-point pressure exceeds the final pressure value, the final flux value is used. This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE.
- fpThe name of the user object for fluid properties
C++ Type:UserObjectName
Controllable:No
Description:The name of the user object for fluid properties
- p_or_t_valsTuple of pressure (or temperature) values. Must be monotonically increasing.
C++ Type:std::vector<Real>
Unit:(no unit assumed)
Controllable:No
Description:Tuple of pressure (or temperature) values. Must be monotonically increasing.
- pressurePressure
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Pressure
- variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on
EnthalpySink
Enthalpy (energy) sink layered on a PorousFlowPolyLineSink: the mass outflow is multiplied by the fluid enthalpy at the local pressure and a constant inlet temperature T_in.
Description
EnthalpySink is a PorousFlowPolyLineSink that represents the energy carried by fluid extracted through a polyline (well) sink. It reuses the base class' pressure-dependent mass outflow rate (a piecewise-linear function of the local porepressure, see p_or_t_vals/fluxes), and multiplies that mass outflow by the specific enthalpy of the fluid, evaluated at the local porepressure and a fixed inlet temperature T_in:
where is the base class' mass outflow rate and is computed by the fp fluid properties user object. EnthalpySink should be applied to the energy (temperature) equation, alongside a companion mass sink (typically a plain PorousFlowPolyLineSink acting on the pressure variable) so that the mass and energy removed from the model are consistent.
Because the outflow is multiplied by a fixed T_in rather than the local temperature, this object is most appropriate for a line sink where the fluid leaves the reservoir at a known/controlled temperature (e.g. a fixed-temperature offtake), as opposed to production at the local reservoir temperature (for which PointEnthalpySink or PorousFlowSquarePulsePointEnthalpySink may be more appropriate).
Example Input Syntax
[DiracKernels<<<{"href": "../../syntax/DiracKernels/index.html"}>>>]
[source_h]
type = EnthalpySink<<<{"description": "Enthalpy (energy) sink layered on a PorousFlowPolyLineSink: the mass outflow is multiplied by the fluid enthalpy at the local pressure and a constant inlet temperature T_in.", "href": "EnthalpySink.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
fp<<<{"description": "The name of the user object for fluid properties"}>>> = simple_fluid
pressure<<<{"description": "Pressure"}>>> = pressure
T_in<<<{"description": "The incoming fluid temperature"}>>> = 320
fluid_phase<<<{"description": "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink. For p_or_t=temperature, and without any use_*, this parameter is irrelevant"}>>> = 0
line_base<<<{"description": "Line base point x,y,z coordinates. This is the same format as a single-line point_file. Note this is only used if there is no point file specified."}>>> = '1 0.5 0.5'
line_direction<<<{"description": "Line direction. Note this is only used if there is only one point in the point_file."}>>> = '0 1 0'
line_length<<<{"description": "Line length. Note this is only used if there is only one point in the point_file."}>>> = 0.3
p_or_t_vals<<<{"description": "Tuple of pressure (or temperature) values. Must be monotonically increasing."}>>> = '0 1e10'
fluxes<<<{"description": "Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed). These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1. A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature). If a quad-point pressure is less than the first pressure value, the first flux value is used. If quad-point pressure exceeds the final pressure value, the final flux value is used. This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE."}>>> = '0.05 0.05'
SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = outflow_energy_uo
[]
[]Input Parameters
- allow_moving_sourcesFalseIf true, allow Dirac sources to move, even if the mesh does not move, during the simulation.
Default:False
C++ Type:bool
Controllable:No
Description:If true, allow Dirac sources to move, even if the mesh does not move, during the simulation.
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- fluid_phase0The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink. For p_or_t=temperature, and without any use_*, this parameter is irrelevant
Default:0
C++ Type:unsigned int
Controllable:No
Description:The fluid phase whose pressure (and potentially mobility, enthalpy, etc) controls the flux to the line sink. For p_or_t=temperature, and without any use_*, this parameter is irrelevant
- function_ofpressureModifying functions will be a function of either pressure and permeability (eg, for boreholes that pump fluids) or temperature and thermal conductivity (eg, for boreholes that pump pure heat with no fluid flow)
Default:pressure
C++ Type:MooseEnum
Options:pressure, temperature
Controllable:No
Description:Modifying functions will be a function of either pressure and permeability (eg, for boreholes that pump fluids) or temperature and thermal conductivity (eg, for boreholes that pump pure heat with no fluid flow)
- line_baseLine base point x,y,z coordinates. This is the same format as a single-line point_file. Note this is only used if there is no point file specified.
C++ Type:std::vector<Real>
Unit:(no unit assumed)
Controllable:No
Description:Line base point x,y,z coordinates. This is the same format as a single-line point_file. Note this is only used if there is no point file specified.
- line_direction0 0 1Line direction. Note this is only used if there is only one point in the point_file.
Default:0 0 1
C++ Type:libMesh::VectorValue<Real>
Unit:(no unit assumed)
Controllable:No
Description:Line direction. Note this is only used if there is only one point in the point_file.
- line_length0Line length. Note this is only used if there is only one point in the point_file.
Default:0
C++ Type:Real
Unit:(no unit assumed)
Range:line_length>=0
Controllable:No
Description:Line length. Note this is only used if there is only one point in the point_file.
- mass_fraction_componentThe index corresponding to a fluid component. If supplied, the flux will be multiplied by the nodal mass fraction for the component
C++ Type:unsigned int
Controllable:No
Description:The index corresponding to a fluid component. If supplied, the flux will be multiplied by the nodal mass fraction for the component
- multiplying_var1.0Fluxes will be moultiplied by this variable
Default:1.0
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Fluxes will be moultiplied by this variable
- point_fileThe file containing the coordinates of the points and their weightings that approximate the line sink. The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii. Each line in the file must contain a space-separated weight and coordinate, viz r x y z. For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure. If your file contains just one point, you must also specify the line_length and line_direction parameters. Note that you will get segementation faults if your points do not lie within your mesh!
C++ Type:std::string
Controllable:No
Description:The file containing the coordinates of the points and their weightings that approximate the line sink. The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii. Each line in the file must contain a space-separated weight and coordinate, viz r x y z. For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure. If your file contains just one point, you must also specify the line_length and line_direction parameters. Note that you will get segementation faults if your points do not lie within your mesh!
- point_not_found_behaviorIGNOREBy default (IGNORE), it is ignored if an added point cannot be located in the specified subdomains. If this option is set to ERROR, this situation will result in an error. If this option is set to WARNING, then a warning will be issued.
Default:IGNORE
C++ Type:MooseEnum
Options:ERROR, WARNING, IGNORE
Controllable:No
Description:By default (IGNORE), it is ignored if an added point cannot be located in the specified subdomains. If this option is set to ERROR, this situation will result in an error. If this option is set to WARNING, then a warning will be issued.
- use_enthalpyFalseMultiply the flux by the fluid enthalpy
Default:False
C++ Type:bool
Controllable:No
Description:Multiply the flux by the fluid enthalpy
- use_internal_energyFalseMultiply the flux by the fluid internal energy
Default:False
C++ Type:bool
Controllable:No
Description:Multiply the flux by the fluid internal energy
- use_mobilityFalseMultiply the flux by the fluid mobility
Default:False
C++ Type:bool
Controllable:No
Description:Multiply the flux by the fluid mobility
- use_relative_permeabilityFalseMultiply the flux by the fluid relative permeability
Default:False
C++ Type:bool
Controllable:No
Description:Multiply the flux by the fluid relative permeability
- weight_reporterreporter weight name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file C++ Type:ReporterName
Controllable:No
Description:reporter weight name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file - x_coord_reporterreporter x-coordinate name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file C++ Type:ReporterName
Controllable:No
Description:reporter x-coordinate name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file - y_coord_reporterreporter y-coordinate name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file C++ Type:ReporterName
Controllable:No
Description:reporter y-coordinate name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file - z_coord_reporterreporter z-coordinate name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file C++ Type:ReporterName
Controllable:No
Description:reporter z-coordinate name of line sink. This uses the reporter syntax
/ . Each point must adhere to the same requirements as those that would be given if using point_file
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
- extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
- matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
- matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill
Contribution To Tagged Field Data Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- drop_duplicate_pointsTrueBy default points added to a DiracKernel are dropped if a point at the same locationhas been added before. If this option is set to false duplicate points are retainedand contribute to residual and Jacobian.
Default:True
C++ Type:bool
Controllable:No
Description:By default points added to a DiracKernel are dropped if a point at the same locationhas been added before. If this option is set to false duplicate points are retainedand contribute to residual and Jacobian.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
Default:nearest_node_connected_sides
C++ Type:MooseEnum
Options:nearest_node_connected_sides, all_proximate_sides
Controllable:No
Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- 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
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
- (test/tests/dirackernels/jacobian/inp.i)
- (test/tests/dirackernels/enthalpy_sink/inp.i)
- (examples/rtes_htm/Cranfield_HTM_push-pull_season_caps_water_newboundary_coarse.i)
- (examples/rtes_benchmark/Cranfield_HT_season_caps_water_oldboundary_noComp_coarse.i)
- (examples/rtes_benchmark/Cranfield_HT_season_caps_water_newboundary_coarse.i)
- (test/tests/dirackernels/jacobian/inp_scaled.i)
- (examples/rtes_htm/Cranfield_HTM_continous_season_caps_water_newboundary_coarse.i)
- (examples/rtes_stochastic/base_cycle_50.i)
- (examples/rtes_benchmark/Cranfield_HT_season_caps_water_oldboundary_coarse.i)
(test/tests/dirackernels/enthalpy_sink/inp.i)
# Analytic check for the EnthalpySink DiracKernel.
#
# EnthalpySink layers an enthalpy multiplier on top of PorousFlowPolyLineSink:
# it removes mass at a rate given by a piecewise-linear function of the local
# porepressure (the usual PorousFlowPolyLineSink behaviour), and multiplies
# that mass outflow by the fluid enthalpy evaluated at (local pressure, T_in).
#
# NOTE on setup: a hand-rolled minimal [Materials] block (mimicking other
# dirackernels tests in this directory) was tried first, but PorousFlowLineSink
# (the base class of EnthalpySink) reports "you do not have a quadpoint
# porepressure material" when EnthalpySink's own required "fp"/"pressure"
# parameters are present, even though the identical [Materials] block works
# fine for a bare PorousFlowPolyLineSink. Building the porous-flow materials
# via the [PorousFlowUnsaturated] action (the pattern used by the working
# examples/rtes_* input files that already use EnthalpySink) avoids the
# problem, so that action is used here instead.
#
# To get an exact hand-derivable result:
# - "fluxes" is a FLAT table ('RATE RATE'), so the mass outflow rate is
# exactly RATE = 0.05 kg/s regardless of the local porepressure.
# - SimpleFluidProperties is used with porepressure_coefficient = 0, which
# forces h_from_p_T(P,T) = e_from_p_T(P,T) = cv*T exactly (see
# SimpleFluidProperties::h_from_p_T), independent of pressure/density.
# - T_in is a fixed constant, so h = cv*T_in = 4000*320 = 1.28e6 J/kg exactly.
#
# Note "fluxes" is a flux PER UNIT LENGTH (kg.m^-1.s^-1), and PorousFlowLineSink
# multiplies it by "line_length" (0.3 m here) to get the actual mass rate in
# kg/s, i.e. mass rate = RATE * line_length = 0.05 * 0.3 = 0.015 kg/s.
#
# Hence the exact energy removed from the domain over one timestep of length
# dt is: RATE * line_length * cv * T_in * dt
# = 0.05 * 0.3 * 4000 * 320 * 100 = 1.92e6 J.
# PorousFlowLineSink (the base class) accumulates exactly this quantity, each
# timestep, into the SumQuantityUO, which PorousFlowPlotQuantity reports as
# "outflow_energy". The "residual" postprocessor below is
# abs(outflow_energy - analytic) and should be ~0 (to solver/roundoff
# tolerance) at every timestep.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pressure]
initial_condition = 1e6
[]
[temperature]
initial_condition = 300
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
density0 = 1000
thermal_expansion = 0
cv = 4000
porepressure_coefficient = 0 # forces h_from_p_T(P,T) = cv*T exactly (see class doc)
[]
[]
[PorousFlowUnsaturated]
porepressure = pressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = simple_fluid
add_darcy_aux = false
[]
[Materials]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.5 0 0 0 2.5 0 0 0 2.5'
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 830.0
density = 2750
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[UserObjects]
[outflow_energy_uo]
type = PorousFlowSumQuantity
[]
[]
[DiracKernels]
[source_h]
type = EnthalpySink
variable = temperature
fp = simple_fluid
pressure = pressure
T_in = 320
fluid_phase = 0
line_base = '1 0.5 0.5'
line_direction = '0 1 0'
line_length = 0.3
p_or_t_vals = '0 1e10'
fluxes = '0.05 0.05'
SumQuantityUO = outflow_energy_uo
[]
[]
[Postprocessors]
[outflow_energy]
type = PorousFlowPlotQuantity
uo = outflow_energy_uo
[]
[residual]
type = ParsedPostprocessor
pp_names = 'outflow_energy'
expression = 'abs(outflow_energy - 0.05 * 0.3 * 4000 * 320 * 100)'
[]
[]
[Preconditioning]
[preferred]
type = SMP
full = true
petsc_options_iname = '-pc_type'
petsc_options_value = ' lu '
[]
[]
[Executioner]
type = Transient
solve_type = Newton
nl_abs_tol = 1e-10
dt = 100
end_time = 300
[]
[Outputs]
csv = true
execute_on = 'initial timestep_end'
[]
(test/tests/dirackernels/jacobian/inp.i)
# Jacobian test for the RTES enthalpy DiracKernels reviewed for missing/incorrect Jacobian terms.
#
# How this is run: the PetscJacobianTester in ./tests injects -snes_test_jacobian itself; it is NOT
# supplied by the [Preconditioning] block below, and -snes_type=test has not existed since PETSc
# 3.9. Because run_sim defaults to False the tester also adds
# -snes_type ksponly -ksp_type preonly -pc_type none -snes_convergence_test skip
# so there is NO nonlinear solve: with Executioner/num_steps=1, exactly one Jacobian is checked, at
# the initial state. One consequence is that only factor = 1.0 of the SquarePulse kernels' six-case
# pulse ramp is ever exercised -- which is why that ramp now lives in a single pulseFactor() helper
# per class rather than being copied into each of computeQpResidual/Jacobian/OffDiagJacobian.
#
# The defects under test:
# - EnthalpySink, FunctionEnthalpySink: missing computeQpBaseOutflowJacobian override (were
# missing both the fluid-enthalpy factor and its dh/dp chain-rule term).
# - PorousFlowSquarePulsePointEnthalpySource: off-diagonal Jacobian was missing a leading minus
# relative to its residual.
# - PointEnthalpySink, PorousFlowSquarePulsePointEnthalpySink, PointEnthalpySourceSinkFromFunction
# (extraction branch): computeQpJacobian() returned 0, discarding a genuine dh/dT diagonal
# term (these are applied to the temperature variable, and h depends on the local
# PorousFlow_temperature_qp).
#
# All of the analytic residual tests elsewhere in test/tests/dirackernels/ deliberately set
# SimpleFluidProperties' porepressure_coefficient=0 to get closed-form residuals -- but that also
# zeroes dh/dp everywhere, and (combined with most of those tests using solve=false) means none of
# them exercise a Jacobian at all. This test leaves porepressure_coefficient at its physical
# default of 1.0 so dh/dp is live.
#
# This input keeps the full PorousFlowUnsaturated physics, which makes it representative but also
# means PorousFlowMatrixInternalEnergy's ~1e5 energy diagonal dominates ||J||_F, leaving too coarse
# an absolute budget to see the O(1e-5) pressure-block entries. The three dh/dT diagonals are
# caught here; the pressure-block terms are covered by the companion inp_scaled.i, which is the
# same eight kernels on an O(1)-scaled problem. See that file's header.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pressure]
[]
[temperature]
[]
[]
[ICs]
[pressure_ic]
type = RandomIC
variable = pressure
min = 1e5
max = 2e6
[]
[temperature_ic]
type = RandomIC
variable = temperature
min = 280
max = 350
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
density0 = 1000
thermal_expansion = 0
cv = 4000
# porepressure_coefficient is left at its physical default (1.0) -- see header comment.
[]
[]
[PorousFlowUnsaturated]
porepressure = pressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = simple_fluid
add_darcy_aux = false
[]
[Materials]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.5 0 0 0 2.5 0 0 0 2.5'
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 830.0
density = 2750
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[UserObjects]
[enthalpy_sink_uo]
type = PorousFlowSumQuantity
[]
[function_enthalpy_sink_uo]
type = PorousFlowSumQuantity
[]
[source_sink_fn_uo]
type = PorousFlowSumQuantity
[]
[mass_source_sink_fn_uo]
type = PorousFlowSumQuantity
[]
[]
[Functions]
[t_in_fn]
type = ParsedFunction
expression = '310 + 5*y'
[]
[mass_flux_fn]
type = ParsedFunction
expression = '0.01'
[]
[temperature_fn]
type = ParsedFunction
expression = '0'
[]
[]
[DiracKernels]
[enthalpy_sink]
type = EnthalpySink
variable = temperature
fp = simple_fluid
pressure = pressure
T_in = 320
fluid_phase = 0
line_base = '1 0.1 0.1 0.1'
line_direction = '0 1 0'
line_length = 0.2
# Deliberately non-flat: a flat table makes _sink_func.sampleDerivative() identically 0, which
# multiplies the h * outflowp half of the chain-rule fix by literal zero.
p_or_t_vals = '0 1e10'
fluxes = '0.04 0.06'
SumQuantityUO = enthalpy_sink_uo
[]
[function_enthalpy_sink]
type = FunctionEnthalpySink
variable = temperature
fp = simple_fluid
pressure = pressure
function = t_in_fn
fluid_phase = 0
line_base = '1 0.6 0.1 0.1'
line_direction = '0 1 0'
line_length = 0.2
p_or_t_vals = '0 1e10'
fluxes = '0.04 0.06'
SumQuantityUO = function_enthalpy_sink_uo
[]
[point_enthalpy_sink]
type = PointEnthalpySink
variable = temperature
mass_flux = mass_flux_pp
fp = simple_fluid
pressure = pressure
point = '0.2 0.8 0.2'
[]
[square_pulse_sink]
type = PorousFlowSquarePulsePointEnthalpySink
variable = temperature
mass_flux = 0.01
fp = simple_fluid
pressure = pressure
point = '0.5 0.8 0.2'
[]
[square_pulse_source]
type = PorousFlowSquarePulsePointEnthalpySource
variable = temperature
mass_flux = 0.01
fp = simple_fluid
pressure = pressure
T_in = 330
point = '0.8 0.8 0.2'
[]
[enthalpy_source_from_function]
# Has a hand-written nonzero pressure off-diagonal and a hard-coded return 0. diagonal; neither
# was verified by any test before this one.
type = PointEnthalpySourceFromFunction
variable = temperature
mass_flux = mass_flux_pp
fp = simple_fluid
pressure = pressure
function = t_in_fn
point = '0.5 0.2 0.8'
[]
[mass_source_sink_from_function]
# Contributes no Jacobian at all (the residual is a pure Function of time and position); this
# pins that down rather than leaving it asserted only in a comment.
type = PointSourceSinkFromFunction
variable = pressure
mass_flux_function = mass_flux_fn
SumQuantityUO = mass_source_sink_fn_uo
point = '0.8 0.2 0.8'
[]
[source_sink_from_function]
type = PointEnthalpySourceSinkFromFunction
variable = temperature
SumQuantityUO = source_sink_fn_uo
mass_flux_function = mass_flux_fn
temperature_function = temperature_fn
pressure = pressure
fp = simple_fluid
point = '0.2 0.2 0.8'
[]
[]
[Postprocessors]
[mass_flux_pp]
type = ConstantPostprocessor
value = 0.01
[]
[]
[Preconditioning]
[check]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 2
[]
[Outputs]
exodus = false
[]
(test/tests/dirackernels/enthalpy_sink/inp.i)
# Analytic check for the EnthalpySink DiracKernel.
#
# EnthalpySink layers an enthalpy multiplier on top of PorousFlowPolyLineSink:
# it removes mass at a rate given by a piecewise-linear function of the local
# porepressure (the usual PorousFlowPolyLineSink behaviour), and multiplies
# that mass outflow by the fluid enthalpy evaluated at (local pressure, T_in).
#
# NOTE on setup: a hand-rolled minimal [Materials] block (mimicking other
# dirackernels tests in this directory) was tried first, but PorousFlowLineSink
# (the base class of EnthalpySink) reports "you do not have a quadpoint
# porepressure material" when EnthalpySink's own required "fp"/"pressure"
# parameters are present, even though the identical [Materials] block works
# fine for a bare PorousFlowPolyLineSink. Building the porous-flow materials
# via the [PorousFlowUnsaturated] action (the pattern used by the working
# examples/rtes_* input files that already use EnthalpySink) avoids the
# problem, so that action is used here instead.
#
# To get an exact hand-derivable result:
# - "fluxes" is a FLAT table ('RATE RATE'), so the mass outflow rate is
# exactly RATE = 0.05 kg/s regardless of the local porepressure.
# - SimpleFluidProperties is used with porepressure_coefficient = 0, which
# forces h_from_p_T(P,T) = e_from_p_T(P,T) = cv*T exactly (see
# SimpleFluidProperties::h_from_p_T), independent of pressure/density.
# - T_in is a fixed constant, so h = cv*T_in = 4000*320 = 1.28e6 J/kg exactly.
#
# Note "fluxes" is a flux PER UNIT LENGTH (kg.m^-1.s^-1), and PorousFlowLineSink
# multiplies it by "line_length" (0.3 m here) to get the actual mass rate in
# kg/s, i.e. mass rate = RATE * line_length = 0.05 * 0.3 = 0.015 kg/s.
#
# Hence the exact energy removed from the domain over one timestep of length
# dt is: RATE * line_length * cv * T_in * dt
# = 0.05 * 0.3 * 4000 * 320 * 100 = 1.92e6 J.
# PorousFlowLineSink (the base class) accumulates exactly this quantity, each
# timestep, into the SumQuantityUO, which PorousFlowPlotQuantity reports as
# "outflow_energy". The "residual" postprocessor below is
# abs(outflow_energy - analytic) and should be ~0 (to solver/roundoff
# tolerance) at every timestep.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pressure]
initial_condition = 1e6
[]
[temperature]
initial_condition = 300
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
density0 = 1000
thermal_expansion = 0
cv = 4000
porepressure_coefficient = 0 # forces h_from_p_T(P,T) = cv*T exactly (see class doc)
[]
[]
[PorousFlowUnsaturated]
porepressure = pressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = simple_fluid
add_darcy_aux = false
[]
[Materials]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.5 0 0 0 2.5 0 0 0 2.5'
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 830.0
density = 2750
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[UserObjects]
[outflow_energy_uo]
type = PorousFlowSumQuantity
[]
[]
[DiracKernels]
[source_h]
type = EnthalpySink
variable = temperature
fp = simple_fluid
pressure = pressure
T_in = 320
fluid_phase = 0
line_base = '1 0.5 0.5'
line_direction = '0 1 0'
line_length = 0.3
p_or_t_vals = '0 1e10'
fluxes = '0.05 0.05'
SumQuantityUO = outflow_energy_uo
[]
[]
[Postprocessors]
[outflow_energy]
type = PorousFlowPlotQuantity
uo = outflow_energy_uo
[]
[residual]
type = ParsedPostprocessor
pp_names = 'outflow_energy'
expression = 'abs(outflow_energy - 0.05 * 0.3 * 4000 * 320 * 100)'
[]
[]
[Preconditioning]
[preferred]
type = SMP
full = true
petsc_options_iname = '-pc_type'
petsc_options_value = ' lu '
[]
[]
[Executioner]
type = Transient
solve_type = Newton
nl_abs_tol = 1e-10
dt = 100
end_time = 300
[]
[Outputs]
csv = true
execute_on = 'initial timestep_end'
[]
(examples/rtes_htm/Cranfield_HTM_push-pull_season_caps_water_newboundary_coarse.i)
#---geothermal battery project - Cranfield - cross validation with TOUGH simulator
#---define initial pressure and temperature
pp_ini_bc = 3e+7 # Pa
T_ini_bc = 399.15 # K
sigmaV_ini_bc = 67.3E6 # Pa
sigmaH_ini_bc = 42.9E6 # Pa
well_length = 1.2 #lower two-thirds of the formation
inj_ext_flux= ${fparse 3/well_length/4 } # 3 kg/s over injection length with 1/4 sysmentrical model
# Darcy flow with heat advection and conduction
[Mesh]
[./fmg]
type = FileMeshGenerator
file = ../rtes_benchmark/Cranfield_caps_new_coarse.e
[]
[]
#############################################################
[GlobalParams]
PorousFlowDictator = dictator
displacements = 'disp_x disp_y disp_z'
[]
####################################################‰########
[PorousFlowUnsaturated]
relative_permeability_type = COREY
relative_permeability_exponent = 0
add_darcy_aux = true
coupling_type = ThermoHydroMechanical
gravity = '0 0 -9.8'
porepressure = porepressure
temperature = temperature
eigenstrain_names = thermal_contribution
fp = tabulated_water
use_displaced_mesh = true
biot_coefficient = 0.9
[]
###########################################################
[Variables]
[./porepressure]
[../]
[./temperature]
scaling = 1E-8
[../]
[./disp_x]
scaling = 1E-10
[../]
[./disp_y]
scaling = 1E-10
[../]
[./disp_z]
scaling = 1E-10
[../]
[]
###########################################################
[AuxVariables]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
initial_condition = -${sigmaH_ini_bc}
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
initial_condition = -${sigmaH_ini_bc}
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
initial_condition = -${sigmaV_ini_bc}
[../]
[./density]
order = CONSTANT
family = MONOMIAL
[../]
[./viscosity]
order = CONSTANT
family = MONOMIAL
[../]
[./porosity]
family = monomial
order = constant
[../]
[./perm_z]
family = monomial
order = constant
[../]
[./perm_x]
family = monomial
order = constant
[../]
[]
##############################################################
[AuxKernels]
[./stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[../]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[../]
[./density]
type = MaterialRealAux
variable = density
property = PorousFlow_fluid_phase_density_qp0
execute_on = TIMESTEP_END
[../]
[./viscosity]
type = MaterialRealAux
variable = viscosity
property = PorousFlow_viscosity_qp0
execute_on = TIMESTEP_END
[../]
[./porosity]
type = MaterialRealAux
property = PorousFlow_porosity_qp
variable = porosity
execute_on = TIMESTEP_END
[../]
[./perm_z]
type = PorousFlowPropertyAux
property = permeability
variable = perm_z
row = 2
column = 2
[../]
[./perm_x]
type = PorousFlowPropertyAux
property = permeability
variable = perm_x
row = 1
column = 1
[../]
[]
############################################################
[UserObjects]
[./fluid_mass_in_inc]
type = PorousFlowSumQuantity
[../]
[./fluid_mass_out_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_in_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_out_inc]
type = PorousFlowSumQuantity
[../]
[]
############################################################
[ICs]
[./init_pp]
type = FunctionIC
function = ${pp_ini_bc}
variable = porepressure
[../]
[./init_temp]
type = ConstantIC
value = ${T_ini_bc}
variable = temperature
[../]
[]
###############################################################
[BCs]
[./roller_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = south
[../]
[./total_stress_y]
type = NeumannBC
variable = disp_y
value = -${sigmaH_ini_bc}
boundary = north
[../]
[./roller_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = west
[../]
[./total_stress_x]
type = NeumannBC
variable = disp_x
value = -${sigmaH_ini_bc}
boundary = east
[../]
[./roller_z]
type = DirichletBC
variable = disp_z
value = 0
boundary = bottom
[../]
[./overburden_total_stress]
type = NeumannBC
variable = disp_z
value = -${sigmaV_ini_bc}
boundary = top
[../]
[]
############################################################
[DiracKernels]
[./summer_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_injection_T]
type = FunctionEnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
function = T_inj_function
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./summer_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_injection_T]
type = FunctionEnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
function = T_inj_function
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./rest_injection_P]
type = PorousFlowPolyLineSink
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = 0
fluxes = 0
[../]
[./rest_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 453.15
fp = tabulated_water
p_or_t_vals = 0
fluxes = 0
[../]
[./rest_production_P]
type = PorousFlowPolyLineSink
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = 0
fluxes = 0
[../]
[./rest_production_T]
type = PorousFlowPolyLineSink
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_enthalpy = true
p_or_t_vals = 0
fluxes = 0
[../]
[]
############################################################
[Controls]
[summer_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_summer
enable_objects = 'DiracKernels::summer_injection_P DiracKernels::summer_injection_T DiracKernels::summer_production_P DiracKernels::summer_production_T'
[]
[winter_extraction]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_winter
enable_objects = 'DiracKernels::winter_injection_P DiracKernels::winter_injection_T DiracKernels::winter_production_P DiracKernels::winter_production_T'
[]
[rest]
type = ConditionalFunctionEnableControl
conditional_function = rest_function
enable_objects = 'DiracKernels::rest_injection_P DiracKernels::rest_injection_T DiracKernels::rest_production_P DiracKernels::rest_production_T'
[]
[]
[Functions]
[./inj_function_summer]
type = ParsedFunction
expression = ' ((t/24/3600/365.25-floor(t/24/3600/365.25))>0.375) & ((t/24/3600/365.25-floor(t/24/3600/365.25)) <=0.625) '
[../]
[./inj_function_winter]
type = ParsedFunction
expression = '((t/24/3600/365.25-floor(t/24/3600/365.25))<=0.125) | ((t/24/3600/365.25-floor(t/24/3600/365.25))>0.875)'
[../]
[./T_inj_function]
type = ParsedFunction
expression = '(453.15+399.15)/2+(453.15-399.15)/2*sin((-0.25+t/24/3600/365.25-floor(t/24/3600/365.25))*2*3.1415926535897932)'
[../]
[./rest_function]
type = ParsedFunction
expression = '((t/24/3600/365.25-floor(t/24/3600/365.25))>0.125 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <= 0.375) | ((t/24/3600/365.25-floor(t/24/3600/365.25))>0.625 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <=0.875) '
[../]
[]
############################################################
[FluidProperties]
[./true_water]
type = Water97FluidProperties
[../]
[./tabulated_water]
type = TabulatedFluidProperties
fp = true_water
allow_fp_and_tabulation = true
temperature_min = 275
temperature_max = 600
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_file = water97_tabulated.csv
[../]
[]
############################################################
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
shear_modulus = 8.38E9
poissons_ratio = 0.225
[../]
[./strain]
type = ComputeSmallStrain
eigenstrain_names = 'thermal_contribution ini_stress'
[../]
[./ini_strain]
type = ComputeEigenstrainFromInitialStress
initial_stress = '-${sigmaH_ini_bc} 0 0 0 -${sigmaH_ini_bc} 0 0 0 -${sigmaV_ini_bc}'
eigenstrain_name = ini_stress
[../]
[./thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
stress_free_temperature = ${T_ini_bc} # using the intial temperature
thermal_expansion_coeff = 1.0e-5 # this is the linear thermal expansion coefficient
eigenstrain_name = thermal_contribution
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./density_aquifer]
type = GenericConstantMaterial
prop_names = density
prop_values = 2600.0
[../]
[./internal_energy_aquifer]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 920.0
density = 2600.0
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./thermal_conductivity_aquifer]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.51 0 0 0 2.51 0 0 0 2.51'
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./porosity_well]
type = PorousFlowPorosityConst
block = 'inj_well ext_well'
porosity = 0.299
[../]
[./permeability_well]
type = PorousFlowPermeabilityConst
block = 'inj_well ext_well'
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-09'
[../]
[./porosity_caps]
type = PorousFlowPorosityConst
block = caps
porosity = 0.01
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-18 0 0 0 1E-18 0 0 0 1E-18'
[../]
[./porosity_layer_1]
type = PorousFlowPorosity
block = layer_1
porosity_zero = 0.143
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc} # using the intial pore presure
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_1]
type = PorousFlowPermeabilityKozenyCarman
block = layer_1
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.3887'
f = 9.9856e-12
d = 1
m = 2
n = 5
[../]
[./porosity_layer_2]
type = PorousFlowPorosity
block = layer_2
porosity_zero = 0.158
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_2]
type = PorousFlowPermeabilityKozenyCarman
block = layer_2
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.3902'
f = 2.4912e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_3]
type = PorousFlowPorosity
block = layer_3
porosity_zero = 0.176
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_3]
type = PorousFlowPermeabilityKozenyCarman
block = layer_3
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.3920'
f = 4.5031e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_4]
type = PorousFlowPorosity
block = layer_4
porosity_zero = 0.188
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_4]
type = PorousFlowPermeabilityKozenyCarman
block = layer_4
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.7292'
f = 5.3904e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_5]
type = PorousFlowPorosity
block = layer_5
porosity_zero = 0.166
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_5]
type = PorousFlowPermeabilityKozenyCarman
block = layer_5
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.7091'
f = 3.0350e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_6]
type = PorousFlowPorosity
block = layer_6
porosity_zero = 0.185
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_6]
type = PorousFlowPermeabilityKozenyCarman
block = layer_6
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.6723'
f = 5.4254e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_7]
type = PorousFlowPorosity
block = layer_7
porosity_zero = 0.24
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_7]
type = PorousFlowPermeabilityKozenyCarman
block = layer_7
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9737'
f = 6.6083e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_8]
type = PorousFlowPorosity
block = layer_8
porosity_zero = 0.26
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_8]
type = PorousFlowPermeabilityKozenyCarman
block = layer_8
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9281'
f = 6.4064e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_9]
type = PorousFlowPorosity
block = layer_9
porosity_zero = 0.285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_9]
type = PorousFlowPermeabilityKozenyCarman
block = layer_9
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9434'
f = 5.7640e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_10]
type = PorousFlowPorosity
block = layer_10
porosity_zero = 0.27
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_10]
type = PorousFlowPermeabilityKozenyCarman
block = layer_10
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9693'
f = 6.0536e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_11]
type = PorousFlowPorosity
block = layer_11
porosity_zero = 0.23
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_11]
type = PorousFlowPermeabilityKozenyCarman
block = layer_11
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.7779'
f = 7.0930e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_12]
type = PorousFlowPorosity
block = layer_12
porosity_zero = 0.306
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_12]
type = PorousFlowPermeabilityKozenyCarman
block = layer_12
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 1'
f = 5.1702e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_13]
type = PorousFlowPorosity
block = layer_13
porosity_zero = .285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_13]
type = PorousFlowPermeabilityKozenyCarman
block = layer_13
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9577'
f = 5.7912e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_14]
type = PorousFlowPorosity
block = layer_14
porosity_zero = .155
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_14]
type = PorousFlowPermeabilityKozenyCarman
block = layer_14
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.0203'
f = 8.6194e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_15]
type = PorousFlowPorosity
block = layer_15
porosity_zero = .299
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_15]
type = PorousFlowPermeabilityKozenyCarman
block = layer_15
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.8759'
f = 5.4697e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_16]
type = PorousFlowPorosity
block = layer_16
porosity_zero = .334
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_16]
type = PorousFlowPermeabilityKozenyCarman
block = layer_16
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9716'
f = 4.5139e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_17]
type = PorousFlowPorosity
block = layer_17
porosity_zero = .32
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_17]
type = PorousFlowPermeabilityKozenyCarman
block = layer_17
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9270'
f = 4.9059e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_18]
type = PorousFlowPorosity
block = layer_18
porosity_zero = .294
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_18]
type = PorousFlowPermeabilityKozenyCarman
block = layer_18
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9752'
f = 5.4915e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_19]
type = PorousFlowPorosity
block = layer_19
porosity_zero = .279
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_19]
type = PorousFlowPermeabilityKozenyCarman
block = layer_19
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.6390'
f = 6.3038e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_20]
type = PorousFlowPorosity
block = layer_20
porosity_zero = .071
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_20]
type = PorousFlowPermeabilityKozenyCarman
block = layer_20
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.0370'
f = 1.2915e-11
d = 1
m = 2
n = 5
[../]
[]
############################################################
[Postprocessors]
[./inlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_in_inc
[../]
[./inlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_in_inc
[../]
[./outlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_out_inc
[../]
[./outlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_out_inc
[../]
[./production_temperature]
type = AverageNodalVariableValue
boundary = extraction
variable = temperature
[../]
[./production_porepressure]
type = AverageNodalVariableValue
boundary = extraction
variable = porepressure
[../]
[./injection_temperature]
type = AverageNodalVariableValue
boundary = injection
variable = temperature
[../]
[./injection_porepressure]
type = AverageNodalVariableValue
boundary = injection
variable = porepressure
[../]
[total]
type = MemoryUsage
mem_units = 'bytes'
execute_on = 'INITIAL TIMESTEP_END'
outputs = none # memory usage is not reproducible; excluded from regression CSV
[]
[]
[Functions]
[./Fiss_Function]
type = PiecewiseLinear
x ='0
3944700 11834100 19723500 27612900
35502300 43391700 51281100 59170500
67059900 74949300 82838700 90728100
98617500 106506900 114396300 122285700
130175100 138064500 145953900 153843300
161732700 169622100 177511500 185400900
193290300 201179700 209069100 216958500
224847900 232737300 240626700 248516100
256405500 264294900 272184300 280073700
287963100 295852500 303741900 311631300
315576000'
y = '0
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
5'
[../]
[]
############################################################
[Preconditioning]
#active = 'ilu_may_use_less_mem'
active = 'superlu'
[./ilu_may_use_less_mem]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels'
petsc_options_value = 'gmres asm ilu NONZERO 2'
[../]
[./original]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm lu NONZERO'
[../]
[./superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'gmres lu superlu_dist'
[../]
[./sub_pc_superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_mat_solver_package'
petsc_options_value = 'gmres asm lu NONZERO superlu_dist'
[../]
[./preferred]
#per Andy
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[../]
[./basic]
type = SMP
full = true
petsc_options = '-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'
[../]
[]
##########################################################
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 315576000
dtmax = 432000 # 5 days
dtmin = 10
[./TimeStepper]
type = IterationAdaptiveDT
dt = 1000.0
optimal_iterations = 15
force_step_every_function_point = true
timestep_limiting_function = Fiss_Function
[../]
line_search = none
# controls for linear iterations
l_max_its = 50
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 50
nl_abs_tol = 1e-6
[]
###########################################################
[Outputs]
[./exodus]
type = Exodus
elemental_as_nodal = true
sync_only = true
sync_times = '0
3944700 11834100 19723500 27612900
35502300 43391700 51281100 59170500
67059900 74949300 82838700 90728100
98617500 106506900 114396300 122285700
130175100 138064500 145953900 153843300
161732700 169622100 177511500 185400900
193290300 201179700 209069100 216958500
224847900 232737300 240626700 248516100
256405500 264294900 272184300 280073700
287963100 295852500 303741900 311631300
315576000'
[../]
[./Console]
type = Console
#perf_log = false
output_linear = false
output_nonlinear = true
time_step_interval = 1
[../]
[./CSV]
type = CSV
time_step_interval = 1
[../]
[]
###########################################################
(examples/rtes_benchmark/Cranfield_HT_season_caps_water_oldboundary_noComp_coarse.i)
#---geothermal battery project - Cranfield - cross validation with TOUGH simulator
#---define initial pressure and temperature
pp_ini_bc = 3e+7 # Pa
T_ini_bc = 399.15 # K
well_length = 15.6 #lower two-thirds of the formation
inj_ext_flux= ${fparse 3/well_length/4 } # 3 kg/s over injection length with 1/4 sysmentrical model
# Darcy flow with heat advection and conduction
[Mesh]
# Cranfield_caps_old_coarse.e no longer exists; remapped onto the current
# reference mesh (see Cranfield_caps_new_coarse.e block/nodeset names below).
[./fmg]
type = FileMeshGenerator
file = Cranfield_caps_new_coarse.e
[]
[]
#############################################################
[GlobalParams]
PorousFlowDictator = dictator
[]
####################################################‰########
[PorousFlowUnsaturated]
relative_permeability_type = COREY
relative_permeability_exponent = 0
add_darcy_aux = true
coupling_type = ThermoHydro
gravity = '0 0 -9.8'
porepressure = porepressure
temperature = temperature
fp = tabulated_water
use_displaced_mesh = false
[]
###########################################################
[Variables]
[./porepressure]
[../]
[./temperature]
scaling = 1E-8
[../]
[]
###########################################################
[AuxVariables]
[./density]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_z]
family = monomial
order = constant
[../]
[./perm_x]
family = monomial
order = constant
[../]
[]
##############################################################
[AuxKernels]
[./density]
type = MaterialRealAux
variable = density
property = PorousFlow_fluid_phase_density_qp0
execute_on = TIMESTEP_END
[../]
[./perm_z]
type = PorousFlowPropertyAux
property = permeability
variable = perm_z
row = 2
column = 2
[../]
[./perm_x]
type = PorousFlowPropertyAux
property = permeability
variable = perm_x
row = 1
column = 1
[../]
[]
############################################################
[UserObjects]
[./fluid_mass_in_inc]
type = PorousFlowSumQuantity
[../]
[./fluid_mass_out_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_in_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_out_inc]
type = PorousFlowSumQuantity
[../]
[]
############################################################
[ICs]
[./init_pp]
type = FunctionIC
function = ${pp_ini_bc}
variable = porepressure
[../]
[./init_temp]
type = ConstantIC
value = ${T_ini_bc}
variable = temperature
[../]
[]
###############################################################
#[BCs]
# [./P_drained]
# type = DirichletBC
# variable = porepressure
# boundary = 'top bottom'
# value = ${pp_ini_bc}
# [../]
# [./T_cont]
# type = DirichletBC
# variable = temperature
# boundary = 'top bottom'
# value = ${T_ini_bc}
# [../]
#[]
############################################################
[DiracKernels]
[./summer_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 453.15
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./summer_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 399.15
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[]
############################################################
[Controls]
[summer_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_summer
enable_objects = 'DiracKernels::summer_injection_P DiracKernels::summer_injection_T DiracKernels::summer_production_P DiracKernels::summer_production_T'
[]
[winter_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_winter
enable_objects = 'DiracKernels::winter_injection_P DiracKernels::winter_injection_T DiracKernels::winter_production_P DiracKernels::winter_production_T'
[]
[]
[Functions]
[./inj_function_summer]
type = ParsedFunction
expression = '(t/24/3600/365-floor(t/24/3600/365))<=0.25'
[../]
[./inj_function_winter]
type = ParsedFunction
expression = '(t/24/3600/365-floor(t/24/3600/365))>0.5 & (t/24/3600/365-floor(t/24/3600/365)) <=0.75'
[../]
[]
############################################################
[FluidProperties]
[./true_water]
type = Water97FluidProperties
[../]
[./tabulated_water]
type = TabulatedFluidProperties
fp = true_water
allow_fp_and_tabulation = true
temperature_min = 275
temperature_max = 600
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_file = water97_tabulated.csv
[../]
[]
############################################################
[Materials]
[./internal_energy_aquifer]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 920.0
density = 2600.0
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./thermal_conductivity_aquifer]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.51 0 0 0 2.51 0 0 0 2.51'
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
# inj_well/ext_well are separate blocks only in the new mesh (the old mesh had
# no dedicated well blocks, since its DiracKernel wells spanned the whole
# formation thickness). Borrowed verbatim from the newboundary input so every
# block in the (new) mesh has a defined porosity/permeability.
[./porosity_well]
type = PorousFlowPorosityConst
block = 'inj_well ext_well'
porosity = 0.299
[../]
[./permeability_well]
type = PorousFlowPermeabilityConst
block = 'inj_well ext_well'
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-09'
[../]
[./porosity_caps]
type = PorousFlowPorosityConst
block = caps
porosity = 0.1
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1.0E-18 0 0 0 1.0E-18 0 0 0 1.0E-18'
[../]
[./porosity_layer_1]
type = PorousFlowPorosityConst
block = layer_1
porosity = 0.143
[../]
[./permeability_layer_1]
type = PorousFlowPermeabilityConst
block = layer_1
permeability = '.8130E-15 0 0 0 .8130E-15 0 0 0 .3160E-15'
[../]
[./porosity_layer_2]
type = PorousFlowPorosityConst
block = layer_2
porosity = 0.158
[../]
[./permeability_layer_2]
type = PorousFlowPermeabilityConst
block = layer_2
permeability = '.3460E-14 0 0 0 .3460E-14 0 0 0 .1350E-14'
[../]
[./porosity_layer_3]
type = PorousFlowPorosityConst
block = layer_3
porosity = 0.176
[../]
[./permeability_layer_3]
type = PorousFlowPermeabilityConst
block = layer_3
permeability = '.1120E-13 0 0 0 .1120E-13 0 0 0 .4390E-14'
[../]
[./porosity_layer_4]
type = PorousFlowPorosityConst
block = layer_4
porosity = 0.188
[../]
[./permeability_layer_4]
type = PorousFlowPermeabilityConst
block = layer_4
permeability = '.1920E-13 0 0 0 .1920E-13 0 0 0 .1400E-13'
[../]
[./porosity_layer_5]
type = PorousFlowPorosityConst
block = layer_5
porosity = 0.166
[../]
[./permeability_layer_5]
type = PorousFlowPermeabilityConst
block = layer_5
permeability = '.5500E-14 0 0 0 .5500E-14 0 0 0 .3900E-14'
[../]
[./porosity_layer_6]
type = PorousFlowPorosityConst
block = layer_6
porosity = 0.185
[../]
[./permeability_layer_6]
type = PorousFlowPermeabilityConst
block = layer_6
permeability = '.1770E-13 0 0 0 .1770E-13 0 0 0 .1190E-13'
[../]
[./porosity_layer_7]
type = PorousFlowPorosityConst
block = layer_7
porosity = 0.24
[../]
[./permeability_layer_7]
type = PorousFlowPermeabilityConst
block = layer_7
permeability = '.9110E-13 0 0 0 .9110E-13 0 0 0 .8870E-13'
[../]
[./porosity_layer_8]
type = PorousFlowPorosityConst
block = layer_8
porosity = 0.26
[../]
[./permeability_layer_8]
type = PorousFlowPermeabilityConst
block = layer_8
permeability = '.1390E-12 0 0 0 .1390E-12 0 0 0 .1290E-12'
[../]
[./porosity_layer_9]
type = PorousFlowPorosityConst
block = layer_9
porosity = 0.285
[../]
[./permeability_layer_9]
type = PorousFlowPermeabilityConst
block = layer_9
permeability = '.2120E-12 0 0 0 .2120E-12 0 0 0 .2000E-12'
[../]
[./porosity_layer_10]
type = PorousFlowPorosityConst
block = layer_10
porosity = 0.27
[../]
[./permeability_layer_10]
type = PorousFlowPermeabilityConst
block = layer_10
permeability = '.1630E-12 0 0 0 .1630E-12 0 0 0 .1580E-12'
[../]
[./porosity_layer_11]
type = PorousFlowPorosityConst
block = layer_11
porosity = 0.23
[../]
[./permeability_layer_11]
type = PorousFlowPermeabilityConst
block = layer_11
permeability = '.7700E-13 0 0 0 .7700E-13 0 0 0 .5990E-13'
[../]
[./porosity_layer_12]
type = PorousFlowPorosityConst
block = layer_12
porosity = 0.306
[../]
[./permeability_layer_12]
type = PorousFlowPermeabilityConst
block = layer_12
permeability = '.2880E-12 0 0 0 .2880E-12 0 0 0 .2880E-12'
[../]
[./porosity_layer_13]
type = PorousFlowPorosityConst
block = layer_13
porosity = .285
[../]
[./permeability_layer_13]
type = PorousFlowPermeabilityConst
block = layer_13
permeability = '.2130E-12 0 0 0 .2130E-12 0 0 0 .2040E-12'
[../]
[./porosity_layer_14]
type = PorousFlowPorosityConst
block = layer_14
porosity = .155
[../]
[./permeability_layer_14]
type = PorousFlowPermeabilityConst
block = layer_14
permeability = '.1080E-13 0 0 0 .1080E-13 0 0 0 .2190E-15'
[../]
[./porosity_layer_15]
type = PorousFlowPorosityConst
block = layer_15
porosity = .299
[../]
[./permeability_layer_15]
type = PorousFlowPermeabilityConst
block = layer_15
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-12'
[../]
[./porosity_layer_16]
type = PorousFlowPorosityConst
block = layer_16
porosity = .334
[../]
[./permeability_layer_16]
type = PorousFlowPermeabilityConst
block = layer_16
permeability = '.4230E-12 0 0 0 .4230E-12 0 0 0 .4110E-12'
[../]
[./porosity_layer_17]
type = PorousFlowPorosityConst
block = layer_17
porosity = .32
[../]
[./permeability_layer_17]
type = PorousFlowPermeabilityConst
block = layer_17
permeability = '.3560E-12 0 0 0 .3560E-12 0 0 0 .3300E-12'
[../]
[./porosity_layer_18]
type = PorousFlowPorosityConst
block = layer_18
porosity = .294
[../]
[./permeability_layer_18]
type = PorousFlowPermeabilityConst
block = layer_18
permeability = '.2420E-12 0 0 0 .2420E-12 0 0 0 .2360E-12'
[../]
[./porosity_layer_19]
type = PorousFlowPorosityConst
block = layer_19
porosity = .279
[../]
[./permeability_layer_19]
type = PorousFlowPermeabilityConst
block = layer_19
permeability = '.2050E-12 0 0 0 .2050E-12 0 0 0 .1310E-12'
[../]
[./porosity_layer_20]
type = PorousFlowPorosityConst
block = layer_20
porosity = .071
[../]
[./permeability_layer_20]
type = PorousFlowPermeabilityConst
block = layer_20
permeability = '.2700E-16 0 0 0 .2700E-16 0 0 0 .1000E-17 '
[../]
[]
############################################################
[Postprocessors]
[./is_summer]
type = FunctionValuePostprocessor
function = inj_function_summer
[../]
[./is_winter]
type = FunctionValuePostprocessor
function = inj_function_winter
[../]
[./inlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_in_inc
[../]
[./inlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_in_inc
[../]
[./outlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_out_inc
[../]
[./outlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_out_inc
[../]
[./production_temperature]
type = AverageNodalVariableValue
boundary = extraction
variable = temperature
[../]
[./production_porepressure]
type = AverageNodalVariableValue
boundary = extraction
variable = porepressure
[../]
[./injection_temperature]
type = AverageNodalVariableValue
boundary = injection
variable = temperature
[../]
[./injection_porepressure]
type = AverageNodalVariableValue
boundary = injection
variable = porepressure
[../]
[total]
type = MemoryUsage
mem_units = 'bytes'
execute_on = 'INITIAL TIMESTEP_END'
outputs = none # memory usage is not reproducible; excluded from regression CSV
[]
[]
############################################################
[Preconditioning]
active = 'ilu_may_use_less_mem'
#active = 'superlu'
[./ilu_may_use_less_mem]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels'
petsc_options_value = 'gmres asm ilu NONZERO 2'
[../]
[./original]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm lu NONZERO'
[../]
[./superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'gmres lu superlu_dist'
[../]
[./sub_pc_superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_mat_solver_package'
petsc_options_value = 'gmres asm lu NONZERO superlu_dist'
[../]
[./preferred]
#per Andy
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[../]
[./basic]
type = SMP
full = true
petsc_options = '-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'
[../]
[]
##########################################################
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 31635000
# steady_state_detection = true
# steady_state_start_time = 864000 # 10 days
# steady_state_tolerance = 1e-6
dtmax = 216000 # 2 days
dtmin = 10
[./TimeStepper]
type = IterationAdaptiveDT
dt = 1000.0
optimal_iterations = 10
[../]
# controls for linear iterations
l_max_its = 100
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 50
nl_abs_tol = 1e-5
[]
###########################################################
[Outputs]
[./exodus]
type = Exodus
# file_base = sequence_out_1
# elemental_as_nodal = true
# overwrite = true
# sequence = true
time_step_interval = 20
[../]
[./Console]
type = Console
#perf_log = false
#output_linear = false
output_nonlinear = true
time_step_interval = 1
[../]
[./CSV]
type = CSV
time_step_interval = 1
[../]
[]
###########################################################
(examples/rtes_benchmark/Cranfield_HT_season_caps_water_newboundary_coarse.i)
#---geothermal battery project - Cranfield - cross validation with TOUGH simulator
#---define initial pressure and temperature
pp_ini_bc = 3e+7 # Pa
T_ini_bc = 399.15 # K
well_length = 1.2 #lower two-thirds of the formation
inj_ext_flux= ${fparse 3/well_length/4 } # 3 kg/s over injection length with 1/4 sysmentrical model
# Darcy flow with heat advection and conduction
[Mesh]
[./fmg]
type = FileMeshGenerator
file = Cranfield_caps_new_coarse.e
[]
[]
#############################################################
[GlobalParams]
PorousFlowDictator = dictator
[]
####################################################‰########
[PorousFlowUnsaturated]
relative_permeability_type = COREY
relative_permeability_exponent = 0
add_darcy_aux = true
coupling_type = ThermoHydro
gravity = '0 0 -9.8'
porepressure = porepressure
temperature = temperature
fp = tabulated_water
use_displaced_mesh = false
[]
###########################################################
[Variables]
[./porepressure]
[../]
[./temperature]
scaling = 1E-8
[../]
[]
###########################################################
[AuxVariables]
[./density]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_z]
family = monomial
order = constant
[../]
[./perm_x]
family = monomial
order = constant
[../]
[]
##############################################################
[AuxKernels]
[./density]
type = MaterialRealAux
variable = density
property = PorousFlow_fluid_phase_density_qp0
execute_on = TIMESTEP_END
[../]
[./perm_z]
type = PorousFlowPropertyAux
property = permeability
variable = perm_z
row = 2
column = 2
[../]
[./perm_x]
type = PorousFlowPropertyAux
property = permeability
variable = perm_x
row = 1
column = 1
[../]
[]
############################################################
[UserObjects]
[./fluid_mass_in_inc]
type = PorousFlowSumQuantity
[../]
[./fluid_mass_out_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_in_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_out_inc]
type = PorousFlowSumQuantity
[../]
[]
############################################################
[ICs]
[./init_pp]
type = FunctionIC
function = ${pp_ini_bc}
variable = porepressure
[../]
[./init_temp]
type = ConstantIC
value = ${T_ini_bc}
variable = temperature
[../]
[]
###############################################################
#[BCs]
# [./P_drained]
# type = DirichletBC
# variable = porepressure
# boundary = 'top bottom'
# value = ${pp_ini_bc}
# [../]
# [./T_cont]
# type = DirichletBC
# variable = temperature
# boundary = 'top bottom'
# value = ${T_ini_bc}
# [../]
#[]
############################################################
[DiracKernels]
[./summer_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 453.15
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./summer_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 399.15
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[]
############################################################
[Controls]
[summer_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_summer
enable_objects = 'DiracKernels::summer_injection_P DiracKernels::summer_injection_T DiracKernels::summer_production_P DiracKernels::summer_production_T'
[]
[winter_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_winter
enable_objects = 'DiracKernels::winter_injection_P DiracKernels::winter_injection_T DiracKernels::winter_production_P DiracKernels::winter_production_T'
[]
[]
[Functions]
[./inj_function_summer]
type = ParsedFunction
expression = '(t/24/3600/365-floor(t/24/3600/365))<=0.25'
[../]
[./inj_function_winter]
type = ParsedFunction
expression = '(t/24/3600/365-floor(t/24/3600/365))>0.5 & (t/24/3600/365-floor(t/24/3600/365)) <=0.75'
[../]
[]
############################################################
[FluidProperties]
[./true_water]
type = Water97FluidProperties
[../]
[./tabulated_water]
type = TabulatedFluidProperties
fp = true_water
allow_fp_and_tabulation = true
temperature_min = 275
temperature_max = 600
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_file = water97_tabulated.csv
[../]
[]
############################################################
[Materials]
[./internal_energy_aquifer]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 920.0
density = 2600.0
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./thermal_conductivity_aquifer]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.51 0 0 0 2.51 0 0 0 2.51'
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./porosity_well]
type = PorousFlowPorosityConst
block = 'inj_well ext_well'
porosity = 0.299
[../]
[./permeability_well]
type = PorousFlowPermeabilityConst
block = 'inj_well ext_well'
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-09'
[../]
[./porosity_caps]
type = PorousFlowPorosityConst
block = caps
porosity = 0.01
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-18 0 0 0 1E-18 0 0 0 1E-18'
[../]
[./porosity_layer_1]
type = PorousFlowPorosity
block = layer_1
porosity_zero = 0.143
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_1]
type = PorousFlowPermeabilityConst
block = layer_1
permeability = '.8130E-15 0 0 0 .8130E-15 0 0 0 .3160E-15'
[../]
[./porosity_layer_2]
type = PorousFlowPorosity
block = layer_2
porosity_zero = 0.158
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_2]
type = PorousFlowPermeabilityConst
block = layer_2
permeability = '.3460E-14 0 0 0 .3460E-14 0 0 0 .1350E-14'
[../]
[./porosity_layer_3]
type = PorousFlowPorosity
block = layer_3
porosity_zero = 0.176
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_3]
type = PorousFlowPermeabilityConst
block = layer_3
permeability = '.1120E-13 0 0 0 .1120E-13 0 0 0 .4390E-14'
[../]
[./porosity_layer_4]
type = PorousFlowPorosity
block = layer_4
porosity_zero = 0.188
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_4]
type = PorousFlowPermeabilityConst
block = layer_4
permeability = '.1920E-13 0 0 0 .1920E-13 0 0 0 .1400E-13'
[../]
[./porosity_layer_5]
type = PorousFlowPorosity
block = layer_5
porosity_zero = 0.166
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_5]
type = PorousFlowPermeabilityConst
block = layer_5
permeability = '.5500E-14 0 0 0 .5500E-14 0 0 0 .3900E-14'
[../]
[./porosity_layer_6]
type = PorousFlowPorosity
block = layer_6
porosity_zero = 0.185
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_6]
type = PorousFlowPermeabilityConst
block = layer_6
permeability = '.1770E-13 0 0 0 .1770E-13 0 0 0 .1190E-13'
[../]
[./porosity_layer_7]
type = PorousFlowPorosity
block = layer_7
porosity_zero = 0.24
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_7]
type = PorousFlowPermeabilityConst
block = layer_7
permeability = '.9110E-13 0 0 0 .9110E-13 0 0 0 .8870E-13'
[../]
[./porosity_layer_8]
type = PorousFlowPorosity
block = layer_8
porosity_zero = 0.26
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_8]
type = PorousFlowPermeabilityConst
block = layer_8
permeability = '.1390E-12 0 0 0 .1390E-12 0 0 0 .1290E-12'
[../]
[./porosity_layer_9]
type = PorousFlowPorosity
block = layer_9
porosity_zero = 0.285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_9]
type = PorousFlowPermeabilityConst
block = layer_9
permeability = '.2120E-12 0 0 0 .2120E-12 0 0 0 .2000E-12'
[../]
[./porosity_layer_10]
type = PorousFlowPorosity
block = layer_10
porosity_zero = 0.27
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_10]
type = PorousFlowPermeabilityConst
block = layer_10
permeability = '.1630E-12 0 0 0 .1630E-12 0 0 0 .1580E-12'
[../]
[./porosity_layer_11]
type = PorousFlowPorosity
block = layer_11
porosity_zero = 0.23
[../]
[./permeability_layer_11]
type = PorousFlowPermeabilityConst
block = layer_11
permeability = '.7700E-13 0 0 0 .7700E-13 0 0 0 .5990E-13'
[../]
[./porosity_layer_12]
type = PorousFlowPorosity
block = layer_12
porosity_zero = 0.306
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_12]
type = PorousFlowPermeabilityConst
block = layer_12
permeability = '.2880E-12 0 0 0 .2880E-12 0 0 0 .2880E-12'
[../]
[./porosity_layer_13]
type = PorousFlowPorosity
block = layer_13
porosity_zero = .285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_13]
type = PorousFlowPermeabilityConst
block = layer_13
permeability = '.2130E-12 0 0 0 .2130E-12 0 0 0 .2040E-12'
[../]
[./porosity_layer_14]
type = PorousFlowPorosity
block = layer_14
porosity_zero = .155
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_14]
type = PorousFlowPermeabilityConst
block = layer_14
permeability = '.1080E-13 0 0 0 .1080E-13 0 0 0 .2190E-15'
[../]
[./porosity_layer_15]
type = PorousFlowPorosity
block = layer_15
porosity_zero = .299
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_15]
type = PorousFlowPermeabilityConst
block = layer_15
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-12'
[../]
[./porosity_layer_16]
type = PorousFlowPorosity
block = layer_16
porosity_zero = .334
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_16]
type = PorousFlowPermeabilityConst
block = layer_16
permeability = '.4230E-12 0 0 0 .4230E-12 0 0 0 .4110E-12'
[../]
[./porosity_layer_17]
type = PorousFlowPorosity
block = layer_17
porosity_zero = .32
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_17]
type = PorousFlowPermeabilityConst
block = layer_17
permeability = '.3560E-12 0 0 0 .3560E-12 0 0 0 .3300E-12'
[../]
[./porosity_layer_18]
type = PorousFlowPorosity
block = layer_18
porosity_zero = .294
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_18]
type = PorousFlowPermeabilityConst
block = layer_18
permeability = '.2420E-12 0 0 0 .2420E-12 0 0 0 .2360E-12'
[../]
[./porosity_layer_19]
type = PorousFlowPorosity
block = layer_19
porosity_zero = .279
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_19]
type = PorousFlowPermeabilityConst
block = layer_19
permeability = '.2050E-12 0 0 0 .2050E-12 0 0 0 .1310E-12'
[../]
[./porosity_layer_20]
type = PorousFlowPorosity
block = layer_20
porosity_zero = .071
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_20]
type = PorousFlowPermeabilityConst
block = layer_20
permeability = '.2700E-16 0 0 0 .2700E-16 0 0 0 .1000E-17 '
[../]
[]
############################################################
[Postprocessors]
[./is_summer]
type = FunctionValuePostprocessor
function = inj_function_summer
[../]
[./is_winter]
type = FunctionValuePostprocessor
function = inj_function_winter
[../]
[./inlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_in_inc
[../]
[./inlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_in_inc
[../]
[./outlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_out_inc
[../]
[./outlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_out_inc
[../]
[./production_temperature]
type = AverageNodalVariableValue
boundary = extraction
variable = temperature
[../]
[./production_porepressure]
type = AverageNodalVariableValue
boundary = extraction
variable = porepressure
[../]
[./injection_temperature]
type = AverageNodalVariableValue
boundary = injection
variable = temperature
[../]
[./injection_porepressure]
type = AverageNodalVariableValue
boundary = injection
variable = porepressure
[../]
[total]
type = MemoryUsage
mem_units = 'bytes'
execute_on = 'INITIAL TIMESTEP_END'
outputs = none # memory usage is not reproducible; excluded from regression CSV
[]
[]
############################################################
[Preconditioning]
active = 'ilu_may_use_less_mem'
#active = 'superlu'
[./ilu_may_use_less_mem]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels'
petsc_options_value = 'gmres asm ilu NONZERO 2'
[../]
[./original]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm lu NONZERO'
[../]
[./superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'gmres lu superlu_dist'
[../]
[./sub_pc_superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_mat_solver_package'
petsc_options_value = 'gmres asm lu NONZERO superlu_dist'
[../]
[./preferred]
#per Andy
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[../]
[./basic]
type = SMP
full = true
petsc_options = '-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'
[../]
[]
##########################################################
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 31635000
# steady_state_detection = true
# steady_state_start_time = 864000 # 10 days
# steady_state_tolerance = 1e-6
dtmax = 216000 # 2 days
dtmin = 10
[./TimeStepper]
type = IterationAdaptiveDT
dt = 1000.0
optimal_iterations = 10
[../]
# controls for linear iterations
l_max_its = 100
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 50
nl_abs_tol = 1e-5
[]
###########################################################
[Outputs]
[./exodus]
type = Exodus
# file_base = sequence_out_1
# elemental_as_nodal = true
# overwrite = true
# sequence = true
time_step_interval = 20
[../]
[./Console]
type = Console
#perf_log = false
#output_linear = false
output_nonlinear = true
time_step_interval = 1
[../]
[./CSV]
type = CSV
time_step_interval = 1
[../]
[]
###########################################################
(test/tests/dirackernels/jacobian/inp_scaled.i)
# Companion to inp.i: the same eight RTES enthalpy DiracKernels, but on a deliberately
# non-dimensionalised problem so that EVERY block of the Jacobian is O(1).
#
# Why a second input exists. PetscJacobianTester gates on ratio_tol, i.e.
# ||J - Jfd||_F / ||J||_F. In inp.i the full PorousFlowUnsaturated action is present, and
# PorousFlowMatrixInternalEnergy (density = 2750, specific_heat_capacity = 830) contributes an
# energy time-derivative diagonal of order 1e5 per node. That single term sets ||J||_F, so the
# absolute error budget at ratio_tol = 1e-7 is ~1e-2 -- while the pressure-block entries these
# kernels actually contribute are mass_flux * dh/dp ~ 1e-5. Three to four orders of magnitude
# separate the two, so a wrong (or missing) pressure off-diagonal is invisible there:
#
# * the dh/dT diagonals of PointEnthalpySink, PorousFlowSquarePulsePointEnthalpySink and
# PointEnthalpySourceSinkFromFunction perturb the norm by ~8 -> caught by inp.i
# * the leading minus of PorousFlowSquarePulsePointEnthalpySource -> ~6e-6, invisible there
# * the dh/dp chain rule of EnthalpySink / FunctionEnthalpySink -> ~6e-6, invisible there
#
# This input follows the pattern of MOOSE's own
# modules/porous_flow/test/tests/jacobian/line_sink02.i: it keeps the PorousFlowDictator and the
# handful of materials the DiracKernels genuinely require, but replaces the action's bulk mass and
# energy kernels with a plain TimeDerivative on each variable, and scales the fluid properties,
# mass fluxes and initial conditions to O(1). With no term dominating, ||J||_F is O(0.1) and the
# pressure-block entries are of the same order, so each defect above moves the ratio by O(1).
#
# porepressure_coefficient is left at its physical default of 1: SimpleFluidProperties computes
# h = e + porepressure_coefficient * p / rho, so setting it to 0 (as the analytic residual tests
# elsewhere in test/tests/dirackernels/ do) would zero dh/dp everywhere and silently remove the
# very term under test.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pressure]
[]
[temperature]
[]
[]
[ICs]
# O(1) ranges, as in line_sink02.i -- these are non-dimensional values chosen to condition the
# Jacobian, not physical pressures and temperatures.
[pressure_ic]
type = RandomIC
variable = pressure
min = 0.2
max = 1.0
[]
[temperature_ic]
type = RandomIC
variable = temperature
min = 1.0
max = 2.0
[]
[]
[Kernels]
# Plain time derivatives in place of the PorousFlow action's mass/energy kernels: they keep the
# diagonal O(1) instead of O(1e5), which is the whole point of this input.
[dummy_pressure]
type = TimeDerivative
variable = pressure
[]
[dummy_temperature]
type = TimeDerivative
variable = temperature
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1.5
density0 = 1
thermal_expansion = 0
viscosity = 1
cv = 1.1
# porepressure_coefficient left at its default of 1 -- see header comment.
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pressure temperature'
number_fluid_phases = 1
number_fluid_components = 1
[]
[enthalpy_sink_uo]
type = PorousFlowSumQuantity
[]
[function_enthalpy_sink_uo]
type = PorousFlowSumQuantity
[]
[source_sink_fn_uo]
type = PorousFlowSumQuantity
[]
[mass_source_sink_fn_uo]
type = PorousFlowSumQuantity
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temperature
[]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = pressure
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1 0 0 0 1 0 0 0 1'
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[Functions]
[t_in_fn]
type = ParsedFunction
expression = '1.4 + 0.2*y'
[]
[mass_flux_fn]
type = ParsedFunction
expression = '0.1'
[]
[temperature_fn]
type = ParsedFunction
expression = '0'
[]
[]
[Postprocessors]
[mass_flux_pp]
type = ConstantPostprocessor
value = 0.1
[]
[]
[DiracKernels]
[enthalpy_sink]
type = EnthalpySink
variable = temperature
fp = simple_fluid
pressure = pressure
T_in = 1.5
fluid_phase = 0
line_base = '1 0.1 0.1 0.1'
line_direction = '0 1 0'
line_length = 0.2
# Deliberately NON-flat: PorousFlowPolyLineSink builds its outflow derivative from
# _sink_func.sampleDerivative(), which is identically 0 for a flat table -- and that would
# multiply the h * outflowp half of the chain-rule fix by literal zero.
p_or_t_vals = '0 2'
fluxes = '0.04 0.06'
SumQuantityUO = enthalpy_sink_uo
[]
[function_enthalpy_sink]
type = FunctionEnthalpySink
variable = temperature
fp = simple_fluid
pressure = pressure
function = t_in_fn
fluid_phase = 0
line_base = '1 0.6 0.1 0.1'
line_direction = '0 1 0'
line_length = 0.2
p_or_t_vals = '0 2'
fluxes = '0.04 0.06'
SumQuantityUO = function_enthalpy_sink_uo
[]
[point_enthalpy_sink]
type = PointEnthalpySink
variable = temperature
mass_flux = mass_flux_pp
fp = simple_fluid
pressure = pressure
point = '0.2 0.8 0.2'
[]
[square_pulse_sink]
type = PorousFlowSquarePulsePointEnthalpySink
variable = temperature
mass_flux = 0.1
fp = simple_fluid
pressure = pressure
point = '0.5 0.8 0.2'
[]
[square_pulse_source]
type = PorousFlowSquarePulsePointEnthalpySource
variable = temperature
mass_flux = 0.1
fp = simple_fluid
pressure = pressure
T_in = 1.7
point = '0.8 0.8 0.2'
[]
[source_sink_from_function]
type = PointEnthalpySourceSinkFromFunction
variable = temperature
SumQuantityUO = source_sink_fn_uo
mass_flux_function = mass_flux_fn
temperature_function = temperature_fn
pressure = pressure
fp = simple_fluid
point = '0.2 0.2 0.8'
[]
[enthalpy_source_from_function]
type = PointEnthalpySourceFromFunction
variable = temperature
mass_flux = mass_flux_pp
fp = simple_fluid
pressure = pressure
function = t_in_fn
point = '0.5 0.2 0.8'
[]
[mass_source_sink_from_function]
type = PointSourceSinkFromFunction
variable = pressure
mass_flux_function = mass_flux_fn
SumQuantityUO = mass_source_sink_fn_uo
point = '0.8 0.2 0.8'
[]
[]
[Preconditioning]
[check]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
exodus = false
[]
(examples/rtes_htm/Cranfield_HTM_continous_season_caps_water_newboundary_coarse.i)
#---geothermal battery project - Cranfield - cross validation with TOUGH simulator
#---define initial pressure and temperature
pp_ini_bc = 3e+7 # Pa
T_ini_bc = 399.15 # K
sigmaV_ini_bc = 67.3E6 # Pa
sigmaH_ini_bc = 42.9E6 # Pa
well_length = 1.2 #lower two-thirds of the formation
inj_ext_flux= ${fparse 3/well_length/4 } # 3 kg/s over injection length with 1/4 sysmentrical model
# Darcy flow with heat advection and conduction
[Mesh]
[./fmg]
type = FileMeshGenerator
file = ../rtes_benchmark/Cranfield_caps_new_coarse.e
[]
[]
#############################################################
[GlobalParams]
PorousFlowDictator = dictator
displacements = 'disp_x disp_y disp_z'
[]
####################################################‰########
[PorousFlowUnsaturated]
relative_permeability_type = COREY
relative_permeability_exponent = 0
add_darcy_aux = true
coupling_type = ThermoHydroMechanical
gravity = '0 0 -9.8'
porepressure = porepressure
temperature = temperature
eigenstrain_names = thermal_contribution
fp = tabulated_water
use_displaced_mesh = true
biot_coefficient = 0.9
[]
###########################################################
[Variables]
[./porepressure]
[../]
[./temperature]
scaling = 1E-8
[../]
[./disp_x]
scaling = 1E-10
[../]
[./disp_y]
scaling = 1E-10
[../]
[./disp_z]
scaling = 1E-10
[../]
[]
###########################################################
[AuxVariables]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
initial_condition = -${sigmaH_ini_bc}
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
initial_condition = -${sigmaH_ini_bc}
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
initial_condition = -${sigmaV_ini_bc}
[../]
[./density]
order = CONSTANT
family = MONOMIAL
[../]
[./viscosity]
order = CONSTANT
family = MONOMIAL
[../]
[./porosity]
family = monomial
order = constant
[../]
[./perm_z]
family = monomial
order = constant
[../]
[./perm_x]
family = monomial
order = constant
[../]
[]
##############################################################
[AuxKernels]
[./stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[../]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[../]
[./density]
type = MaterialRealAux
variable = density
property = PorousFlow_fluid_phase_density_qp0
execute_on = TIMESTEP_END
[../]
[./viscosity]
type = MaterialRealAux
variable = viscosity
property = PorousFlow_viscosity_qp0
execute_on = TIMESTEP_END
[../]
[./porosity]
type = MaterialRealAux
property = PorousFlow_porosity_qp
variable = porosity
execute_on = TIMESTEP_END
[../]
[./perm_z]
type = PorousFlowPropertyAux
property = permeability
variable = perm_z
row = 2
column = 2
[../]
[./perm_x]
type = PorousFlowPropertyAux
property = permeability
variable = perm_x
row = 1
column = 1
[../]
[]
############################################################
[UserObjects]
[./fluid_mass_in_inc]
type = PorousFlowSumQuantity
[../]
[./fluid_mass_out_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_in_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_out_inc]
type = PorousFlowSumQuantity
[../]
[]
############################################################
[ICs]
[./init_pp]
type = FunctionIC
function = ${pp_ini_bc}
variable = porepressure
[../]
[./init_temp]
type = ConstantIC
value = ${T_ini_bc}
variable = temperature
[../]
[]
###############################################################
[BCs]
[./roller_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = south
[../]
[./total_stress_y]
type = NeumannBC
variable = disp_y
value = -${sigmaH_ini_bc}
boundary = north
[../]
[./roller_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = west
[../]
[./total_stress_x]
type = NeumannBC
variable = disp_x
value = -${sigmaH_ini_bc}
boundary = east
[../]
[./roller_z]
type = DirichletBC
variable = disp_z
value = 0
boundary = bottom
[../]
[./overburden_total_stress]
type = NeumannBC
variable = disp_z
value = -${sigmaV_ini_bc}
boundary = top
[../]
[]
############################################################
[DiracKernels]
[./summer_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_injection_T]
type = FunctionEnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
function = T_inj_function
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./summer_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./rest_injection_P]
type = PorousFlowPolyLineSink
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = 0
fluxes = 0
[../]
[./rest_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 453.15
fp = tabulated_water
p_or_t_vals = 0
fluxes = 0
[../]
[./rest_production_P]
type = PorousFlowPolyLineSink
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = 0
fluxes = 0
[../]
[./rest_production_T]
type = PorousFlowPolyLineSink
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 3.6'
line_length = ${well_length}
line_direction = '0 0 1'
use_enthalpy = true
p_or_t_vals = 0
fluxes = 0
[../]
[]
############################################################
[Controls]
[summer_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_summer_winter
enable_objects = 'DiracKernels::summer_injection_P DiracKernels::summer_injection_T DiracKernels::summer_production_P DiracKernels::summer_production_T'
[]
[rest]
type = ConditionalFunctionEnableControl
conditional_function = rest_function
enable_objects = 'DiracKernels::rest_injection_P DiracKernels::rest_injection_T DiracKernels::rest_production_P DiracKernels::rest_production_T'
[]
[]
[Functions]
[./inj_function_summer_winter]
type = ParsedFunction
expression = '((t/24/3600/365.25-floor(t/24/3600/365.25))<=0.125) | ((t/24/3600/365.25-floor(t/24/3600/365.25))>0.375 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <=0.625) | ((t/24/3600/365.25-floor(t/24/3600/365.25))>0.875)'
[../]
[./T_inj_function]
type = ParsedFunction
expression = '(453.15+399.15)/2+(453.15-399.15)/2*sin((-0.25+t/24/3600/365.25-floor(t/24/3600/365.25))*2*3.1415926535897932)'
[../]
[./rest_function]
type = ParsedFunction
expression = '((t/24/3600/365.25-floor(t/24/3600/365.25))>0.125 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <= 0.375) | ((t/24/3600/365.25-floor(t/24/3600/365.25))>0.625 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <=0.875) '
[../]
[]
############################################################
[FluidProperties]
[./true_water]
type = Water97FluidProperties
[../]
[./tabulated_water]
type = TabulatedFluidProperties
fp = true_water
allow_fp_and_tabulation = true
temperature_min = 275
temperature_max = 600
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_file = water97_tabulated.csv
[../]
[]
############################################################
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
shear_modulus = 8.38E9
poissons_ratio = 0.225
[../]
[./strain]
type = ComputeSmallStrain
eigenstrain_names = 'thermal_contribution ini_stress'
[../]
[./ini_strain]
type = ComputeEigenstrainFromInitialStress
initial_stress = '-${sigmaH_ini_bc} 0 0 0 -${sigmaH_ini_bc} 0 0 0 -${sigmaV_ini_bc}'
eigenstrain_name = ini_stress
[../]
[./thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
stress_free_temperature = ${T_ini_bc} # using the intial temperature
thermal_expansion_coeff = 1.0e-5 # this is the linear thermal expansion coefficient
eigenstrain_name = thermal_contribution
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./density_aquifer]
type = GenericConstantMaterial
prop_names = density
prop_values = 2600.0
[../]
[./internal_energy_aquifer]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 920.0
density = 2600.0
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./thermal_conductivity_aquifer]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.51 0 0 0 2.51 0 0 0 2.51'
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./porosity_well]
type = PorousFlowPorosityConst
block = 'inj_well ext_well'
porosity = 0.299
[../]
[./permeability_well]
type = PorousFlowPermeabilityConst
block = 'inj_well ext_well'
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-09'
[../]
[./porosity_caps]
type = PorousFlowPorosityConst
block = caps
porosity = 0.01
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-18 0 0 0 1E-18 0 0 0 1E-18'
[../]
[./porosity_layer_1]
type = PorousFlowPorosity
block = layer_1
porosity_zero = 0.143
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc} # using the intial pore presure
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_1]
type = PorousFlowPermeabilityKozenyCarman
block = layer_1
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.3887'
f = 9.9856e-12
d = 1
m = 2
n = 5
[../]
[./porosity_layer_2]
type = PorousFlowPorosity
block = layer_2
porosity_zero = 0.158
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_2]
type = PorousFlowPermeabilityKozenyCarman
block = layer_2
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.3902'
f = 2.4912e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_3]
type = PorousFlowPorosity
block = layer_3
porosity_zero = 0.176
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_3]
type = PorousFlowPermeabilityKozenyCarman
block = layer_3
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.3920'
f = 4.5031e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_4]
type = PorousFlowPorosity
block = layer_4
porosity_zero = 0.188
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_4]
type = PorousFlowPermeabilityKozenyCarman
block = layer_4
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.7292'
f = 5.3904e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_5]
type = PorousFlowPorosity
block = layer_5
porosity_zero = 0.166
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_5]
type = PorousFlowPermeabilityKozenyCarman
block = layer_5
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.7091'
f = 3.0350e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_6]
type = PorousFlowPorosity
block = layer_6
porosity_zero = 0.185
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_6]
type = PorousFlowPermeabilityKozenyCarman
block = layer_6
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.6723'
f = 5.4254e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_7]
type = PorousFlowPorosity
block = layer_7
porosity_zero = 0.24
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_7]
type = PorousFlowPermeabilityKozenyCarman
block = layer_7
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9737'
f = 6.6083e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_8]
type = PorousFlowPorosity
block = layer_8
porosity_zero = 0.26
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_8]
type = PorousFlowPermeabilityKozenyCarman
block = layer_8
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9281'
f = 6.4064e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_9]
type = PorousFlowPorosity
block = layer_9
porosity_zero = 0.285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_9]
type = PorousFlowPermeabilityKozenyCarman
block = layer_9
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9434'
f = 5.7640e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_10]
type = PorousFlowPorosity
block = layer_10
porosity_zero = 0.27
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_10]
type = PorousFlowPermeabilityKozenyCarman
block = layer_10
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9693'
f = 6.0536e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_11]
type = PorousFlowPorosity
block = layer_11
porosity_zero = 0.23
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_11]
type = PorousFlowPermeabilityKozenyCarman
block = layer_11
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.7779'
f = 7.0930e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_12]
type = PorousFlowPorosity
block = layer_12
porosity_zero = 0.306
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_12]
type = PorousFlowPermeabilityKozenyCarman
block = layer_12
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 1'
f = 5.1702e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_13]
type = PorousFlowPorosity
block = layer_13
porosity_zero = .285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_13]
type = PorousFlowPermeabilityKozenyCarman
block = layer_13
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9577'
f = 5.7912e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_14]
type = PorousFlowPorosity
block = layer_14
porosity_zero = .155
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_14]
type = PorousFlowPermeabilityKozenyCarman
block = layer_14
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.0203'
f = 8.6194e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_15]
type = PorousFlowPorosity
block = layer_15
porosity_zero = .299
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_15]
type = PorousFlowPermeabilityKozenyCarman
block = layer_15
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.8759'
f = 5.4697e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_16]
type = PorousFlowPorosity
block = layer_16
porosity_zero = .334
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_16]
type = PorousFlowPermeabilityKozenyCarman
block = layer_16
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9716'
f = 4.5139e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_17]
type = PorousFlowPorosity
block = layer_17
porosity_zero = .32
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_17]
type = PorousFlowPermeabilityKozenyCarman
block = layer_17
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9270'
f = 4.9059e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_18]
type = PorousFlowPorosity
block = layer_18
porosity_zero = .294
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_18]
type = PorousFlowPermeabilityKozenyCarman
block = layer_18
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.9752'
f = 5.4915e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_19]
type = PorousFlowPorosity
block = layer_19
porosity_zero = .279
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_19]
type = PorousFlowPermeabilityKozenyCarman
block = layer_19
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.6390'
f = 6.3038e-11
d = 1
m = 2
n = 5
[../]
[./porosity_layer_20]
type = PorousFlowPorosity
block = layer_20
porosity_zero = .071
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
mechanical = true
thermal_expansion_coeff = 1.0e-5 # Kim and Hosseini 2013
thermal = true
reference_temperature = ${T_ini_bc} # using the intial temperature
[../]
[./permeability_layer_20]
type = PorousFlowPermeabilityKozenyCarman
block = layer_20
poroperm_function = kozeny_carman_fd2
k_anisotropy = '1 0 0 0 1 0 0 0 0.0370'
f = 1.2915e-11
d = 1
m = 2
n = 5
[../]
[]
############################################################
[Postprocessors]
[./inlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_in_inc
[../]
[./inlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_in_inc
[../]
[./outlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_out_inc
[../]
[./outlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_out_inc
[../]
[./production_temperature]
type = AverageNodalVariableValue
boundary = extraction
variable = temperature
[../]
[./production_porepressure]
type = AverageNodalVariableValue
boundary = extraction
variable = porepressure
[../]
[./injection_temperature]
type = AverageNodalVariableValue
boundary = injection
variable = temperature
[../]
[./injection_porepressure]
type = AverageNodalVariableValue
boundary = injection
variable = porepressure
[../]
[total]
type = MemoryUsage
mem_units = 'bytes'
execute_on = 'INITIAL TIMESTEP_END'
outputs = none # memory usage is not reproducible; excluded from regression CSV
[]
[]
[Functions]
[./Fiss_Function]
type = PiecewiseLinear
x ='0
3944700 11834100 19723500 27612900
35502300 43391700 51281100 59170500
67059900 74949300 82838700 90728100
98617500 106506900 114396300 122285700
130175100 138064500 145953900 153843300
161732700 169622100 177511500 185400900
193290300 201179700 209069100 216958500
224847900 232737300 240626700 248516100
256405500 264294900 272184300 280073700
287963100 295852500 303741900 311631300
315576000'
y = '0
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
5'
[../]
[]
############################################################
[Preconditioning]
#active = 'ilu_may_use_less_mem'
active = 'superlu'
[./ilu_may_use_less_mem]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels'
petsc_options_value = 'gmres asm ilu NONZERO 2'
[../]
[./original]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm lu NONZERO'
[../]
[./superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'gmres lu superlu_dist'
[../]
[./sub_pc_superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_mat_solver_package'
petsc_options_value = 'gmres asm lu NONZERO superlu_dist'
[../]
[./preferred]
#per Andy
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[../]
[./basic]
type = SMP
full = true
petsc_options = '-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'
[../]
[]
##########################################################
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 315576000
dtmax = 432000 # 5 days
dtmin = 10
[./TimeStepper]
type = IterationAdaptiveDT
dt = 1000.0
optimal_iterations = 15
force_step_every_function_point = true
timestep_limiting_function = Fiss_Function
[../]
line_search = none
# controls for linear iterations
l_max_its = 50
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 50
nl_abs_tol = 1e-6
[]
###########################################################
[Outputs]
[./exodus]
type = Exodus
elemental_as_nodal = true
sync_only = true
sync_times = '0
3944700 11834100 19723500 27612900
35502300 43391700 51281100 59170500
67059900 74949300 82838700 90728100
98617500 106506900 114396300 122285700
130175100 138064500 145953900 153843300
161732700 169622100 177511500 185400900
193290300 201179700 209069100 216958500
224847900 232737300 240626700 248516100
256405500 264294900 272184300 280073700
287963100 295852500 303741900 311631300
315576000'
[../]
[./Console]
type = Console
#perf_log = false
output_linear = false
output_nonlinear = true
time_step_interval = 1
[../]
[./CSV]
type = CSV
time_step_interval = 1
[../]
[]
###########################################################
(examples/rtes_stochastic/base_cycle_50.i)
#---geothermal battery project - looking for sweet spot
#---Define two scalar variables to sweep consistently
depth = 2000
pp_ini_bc = ${fparse if(depth<=2000, 1000*9.8*depth, 1000*9.8*depth+(depth-2000)*9.8*1000)}
T_ini_bc = ${fparse (80+3.28084*depth*1.1/100-32)*5/9 + 273.15}
distance_between_wells = 200
pro_well_x = ${fparse 0.1 + distance_between_wells / 2}
inj_well_x = ${fparse 0.1 - distance_between_wells / 2}
well_length = 50 #fixed number in accordance with the mesh
aquifer_mid= ${fparse 0 + 0}
aquifer_top= ${fparse 0 + well_length/2}
aquifer_bot= ${fparse 0 - well_length/2}
#---Define a scalar variable to sweep
inj_ext_exponent = 1.5
inj_ext_flux= ${fparse 10^inj_ext_exponent / well_length / 2 }
injection_temp = 473.15
#---Define a scalar variable to replace the tensor components
Tcond_aquifer = 2.0
#---Define a scalar variable to replace the tensor components
perm_exponent = -12
perm_aquifer = ${fparse 10^perm_exponent}
[Controls]
[summer_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_summer
enable_objects = 'DiracKernel::charge_injection_P DiracKernel::charge_injection_T DiracKernel::charge_production_P DiracKernel::charge_production_T'
[]
[winter_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_winter
enable_objects = 'DiracKernel::recover_injection_P DiracKernel::recover_injection_T DiracKernel::recover_production_P DiracKernel::recover_production_T'
[]
[rest_injection]
type = ConditionalFunctionEnableControl
conditional_function = rest_function
enable_objects = 'DiracKernel::rest_injection_P DiracKernel::rest_injection_T DiracKernel::rest_production_P DiracKernel::rest_production_T'
[]
[]
[Functions]
[./inj_function_summer]
type = ParsedFunction
expression = '(t/24/3600/365.25-floor(t/24/3600/365.25))>0 & (t/24/3600/365.25-floor(t/24/3600/365.25))<=0.25'
[../]
[./inj_function_winter]
type = ParsedFunction
symbol_values = 'termination'
symbol_names = 'b_switch'
expression = '((t/24/3600/365.25-floor(t/24/3600/365.25))>0.5 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <=0.75) & b_switch <= 0'
[../]
[./rest_function]
type = ParsedFunction
symbol_values = 'termination'
symbol_names = 'b_switch'
expression = '(((t/24/3600/365.25-floor(t/24/3600/365.25))>0.25 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <=0.5) | ((t/24/3600/365.25-floor(t/24/3600/365.25))>0.75 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <=1)) | (((t/24/3600/365.25-floor(t/24/3600/365.25))>0.5 & (t/24/3600/365.25-floor(t/24/3600/365.25)) <=0.75) & b_switch >=1)'
[../]
[./Fiss_Function]
type = PiecewiseLinear
x ='0. 7889400. 15778800. 23668200.
31557600. 39447000. 47336400. 55225800.
63115200. 71004600. 78894000. 86783400.
94672800. 102562200. 110451600. 118341000.
126230400. 134119800. 142009200. 149898600.
157788000. 165677400. 173566800. 181456200.
189345600. 197235000. 205124400. 213013800.
220903200. 228792600. 236682000. 244571400.
252460800. 260350200. 268239600. 276129000.
284018400. 291907800. 299797200. 307686600.
315576000.'
y = '0 1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4'
[../]
[]
# Darcy flow with heat advection and conduction
[Mesh]
# The original base_final_50.e had one or more element-block IDs that overflowed the
# unsigned-short block-ID field the current libMesh ExodusII_IO reader uses
# ("restrict_int failed: <id> does not fit in type t"). base_final_50_blockid_fix.e is that
# same mesh with its element-block IDs renumbered to fit (geometry, coordinates, connectivity,
# and block/nodeset/sideset NAMES -- which is what this input and its Materials blocks
# reference by name, not by numeric ID -- are unchanged by a renumbering).
#
# The original base_final_50.e is not present in this repository (on this branch or on
# devel), so the renumbering is not independently auditable from git history alone. The
# renumbered file's own embedded NetCDF metadata is the audit trail instead: its `title`
# global attribute reads
# cubit(s/jinw-mac/projects/falcon/Applied_Energy/base_final_50.e): 10/02/2020: 22
# (i.e. it was produced by Cubit from a file of that name), and its current eb_prop1/eb_names
# (`ncdump -v eb_prop1,eb_names base_final_50_blockid_fix.e`) are:
# caps_HEX8=1, caps_WEDGE=3, aquifer_HEX8=2, aquifer_WEDGE=4
# all comfortably within the unsigned-short range, consistent with a renumbering rather than a
# from-scratch remesh.
[./fmg]
type = FileMeshGenerator
file = base_final_50_blockid_fix.e
[]
[]
#############################################################
[GlobalParams]
PorousFlowDictator = dictator
[]
############################################################
[PorousFlowUnsaturated]
relative_permeability_type = COREY
relative_permeability_exponent = 0
add_darcy_aux = true
coupling_type = ThermoHydro
gravity = '0 0 -9.8'
porepressure = porepressure
temperature = temperature
fp = tabulated_water
use_displaced_mesh = false
[]
###########################################################
[Variables]
[./porepressure]
[../]
[./temperature]
scaling = 1E-6
[../]
[]
###########################################################
[AuxVariables]
[./density]
order = CONSTANT
family = MONOMIAL
[../]
[]
##############################################################
[AuxKernels]
[./density]
type = MaterialRealAux
variable = density
property = PorousFlow_fluid_phase_density_qp0
execute_on = TIMESTEP_END
[../]
[]
############################################################
[UserObjects]
[./fluid_mass_in_inc]
type = PorousFlowSumQuantity
[../]
[./fluid_mass_out_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_in_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_out_inc]
type = PorousFlowSumQuantity
[../]
[]
############################################################
[ICs]
[./init_pp]
type = FunctionIC
function = ${pp_ini_bc}
variable = porepressure
[../]
[./init_temp]
type = ConstantIC
value = ${T_ini_bc}
variable = temperature
[../]
[]
###############################################################
[BCs]
[./P_drained]
type = DirichletBC
variable = porepressure
boundary = 'east west north'
value = ${pp_ini_bc}
[../]
[./T_cont]
type = DirichletBC
variable = temperature
boundary = 'east west north'
value = ${T_ini_bc}
[../]
[]
############################################################
[DiracKernels]
[./charge_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 ${inj_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}' # ~5 kg/s over length of 10(injection_length)/2
[../]
[./charge_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 ${inj_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = ${injection_temp}
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./charge_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 ${pro_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./charge_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 ${pro_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./recover_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 ${inj_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}' # ~5 kg/s over length of 10(injection_length)/2
[../]
[./recover_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 ${inj_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./recover_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 ${pro_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./recover_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 ${pro_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = ${T_ini_bc}
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./rest_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 ${inj_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = '-1e9 1e9'
fluxes = '0 0' # ~5 kg/s over length of 10(injection_length)/2
[../]
[./rest_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 ${inj_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '0 0'
[../]
[./rest_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 ${pro_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
p_or_t_vals = '-1e9 1e9'
fluxes = '0 0'
[../]
[./rest_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 ${pro_well_x} 0 ${aquifer_bot}'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = ${T_ini_bc}
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '0 0'
[../]
[]
############################################################
[FluidProperties]
[./true_water]
type = Water97FluidProperties
[../]
[./tabulated_water]
type = TabulatedFluidProperties
fp = true_water
allow_fp_and_tabulation = true
temperature_min = 275
temperature_max = 600
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_file = water97_tabulated.csv
[../]
[]
############################################################
[Materials]
[./internal_energy_aquifer]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 930.0
density = 2650.0
block = 'aquifer_HEX8 aquifer_WEDGE'
[../]
[./thermal_conductivity_aquifer]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '${Tcond_aquifer} 0 0 0 ${Tcond_aquifer} 0 0 0 ${Tcond_aquifer}'
block = 'aquifer_HEX8 aquifer_WEDGE'
[../]
[./porosity_aquifer]
type = PorousFlowPorosityConst
block = 'aquifer_HEX8 aquifer_WEDGE'
porosity = 0.01
[../]
[./permeability_aquifer]
type = PorousFlowPermeabilityConst
block = 'aquifer_HEX8 aquifer_WEDGE'
permeability = '${perm_aquifer} 0 0 0 ${perm_aquifer} 0 0 0 ${perm_aquifer}'
[../]
[./internal_energy_caps]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 1000.0
density = 2500.0
block = 'caps_HEX8 caps_WEDGE'
[../]
[./thermal_conductivity_caps]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.5 0 0 0 2.5 0 0 0 2.5'
block = 'caps_HEX8 caps_WEDGE'
[../]
[./porosity_caps]
type = PorousFlowPorosityConst
block = 'caps_HEX8 caps_WEDGE'
porosity = 0.01
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = 'caps_HEX8 caps_WEDGE'
permeability = '1E-18 0 0 0 1E-18 0 0 0 1E-18'
[../]
[]
############################################################
[Preconditioning]
active = 'ilu_may_use_less_mem'
#active = 'superlu'
[./ilu_may_use_less_mem]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels'
petsc_options_value = 'gmres asm ilu NONZERO 2'
[../]
[./original]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm lu NONZERO'
[../]
[./superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'gmres lu superlu_dist'
[../]
[./sub_pc_superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_mat_solver_package'
petsc_options_value = 'gmres asm lu NONZERO superlu_dist'
[../]
[./preferred]
#per Andy
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[../]
[./basic]
type = SMP
full = true
petsc_options = '-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'
[../]
[]
############################################################
[Postprocessors]
[./hotwell_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_in_inc
[../]
[./hotwell_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_in_inc
[../]
[./coldwell_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_out_inc
[../]
[./coldwell_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_out_inc
[../]
[./time]
type = TimePostprocessor
execute_on = 'timestep_end'
[]
[./step_dt]
type = TimestepSize
[../]
[./is_winter]
type = FunctionValuePostprocessor
function = inj_function_winter
[../]
[./is_summer]
type = FunctionValuePostprocessor
function = inj_function_summer
[../]
[./total_recovery_time]
type = PorousFlowRecoveryTimeAccumulator
targetpostprocessor = is_winter
dtpostprocessor = step_dt
[../]
[./termination]
type = PorousFlowTemperatureDropTerminator
enthalpypostprocessor = hotwell_enthalpy_J
masspostprocessor = hotwell_mass_kg
T_inj = ${injection_temp}
T_init = ${T_ini_bc}
P_drop = 80
timepostprocessor = time
[../]
[./recovery_rate]
type = PorousFlowRecoveryRateSeason
hotwellenergy = hotwell_enthalpy_J
coldwellenergy = coldwell_enthalpy_J
InjectionIndicator = is_summer
ProductionIndicator = is_winter
[../]
[./total_energy]
type = PorousFlowEnergyAccumulator
hotwellenergy = hotwell_enthalpy_J
coldwellenergy = coldwell_enthalpy_J
ProductionIndicator = is_winter
[../]
[./inj_P_mid]
type = PointValue
execute_on = 'initial timestep_end'
point = '${inj_well_x} 0 ${aquifer_mid}'
variable = porepressure
[../]
[./inj_T_mid]
type = PointValue
execute_on = 'initial timestep_end'
point = '${inj_well_x} 0 ${aquifer_mid}'
variable = temperature
[../]
[./pro_P_mid]
type = PointValue
execute_on = 'initial timestep_end'
point = '${pro_well_x} 0 ${aquifer_mid}'
variable = porepressure
[../]
[./pro_T_mid]
type = PointValue
execute_on = 'initial timestep_end'
point = '${pro_well_x} 0 ${aquifer_mid}'
variable = temperature
[../]
[]
##########################################################
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 315576000
# steady_state_detection = true
# steady_state_start_time = 864000 # 10 days
# steady_state_tolerance = 1e-6
dtmax = 864000 # 10 days
dtmin = 100
[./TimeStepper]
type = IterationAdaptiveDT
dt = 1000.0
optimal_iterations = 10
force_step_every_function_point = true
timestep_limiting_function = Fiss_Function
# timestep_limiting_postprocessor = timestep_pp
[../]
# controls for linear iterations
l_max_its = 50
l_tol = 1e-3
# controls for nonlinear iterations
nl_max_its = 20
nl_rel_tol = 1e-5
nl_abs_tol = 1e-4
[]
###########################################################
[Outputs]
#[./exodus]
# type = Exodus
# file_base = sequence_out_1
# elemental_as_nodal = true
# overwrite = true
# sequence = true
# time_step_interval = 2
#[../]
#[./Console]
# type = Console
# #perf_log = false
# #output_linear = true
# output_nonlinear = true
# time_step_interval = 1
#[../]
[./CSV]
type = CSV
time_step_interval = 1
[../]
[]
###########################################################
(examples/rtes_benchmark/Cranfield_HT_season_caps_water_oldboundary_coarse.i)
#---geothermal battery project - Cranfield - cross validation with TOUGH simulator
#---define initial pressure and temperature
pp_ini_bc = 3e+7 # Pa
T_ini_bc = 399.15 # K
well_length = 15.6 #lower two-thirds of the formation
inj_ext_flux= ${fparse 3/well_length/4 } # 3 kg/s over injection length with 1/4 sysmentrical model
# Darcy flow with heat advection and conduction
[Mesh]
# Cranfield_caps_old_coarse.e no longer exists; remapped onto the current
# reference mesh (see Cranfield_caps_new_coarse.e block/nodeset names below).
[./fmg]
type = FileMeshGenerator
file = Cranfield_caps_new_coarse.e
[]
[]
#############################################################
[GlobalParams]
PorousFlowDictator = dictator
[]
####################################################‰########
[PorousFlowUnsaturated]
relative_permeability_type = COREY
relative_permeability_exponent = 0
add_darcy_aux = true
coupling_type = ThermoHydro
gravity = '0 0 -9.8'
porepressure = porepressure
temperature = temperature
fp = tabulated_water
use_displaced_mesh = false
[]
###########################################################
[Variables]
[./porepressure]
[../]
[./temperature]
scaling = 1E-8
[../]
[]
###########################################################
[AuxVariables]
[./density]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_z]
family = monomial
order = constant
[../]
[./perm_x]
family = monomial
order = constant
[../]
[]
##############################################################
[AuxKernels]
[./density]
type = MaterialRealAux
variable = density
property = PorousFlow_fluid_phase_density_qp0
execute_on = TIMESTEP_END
[../]
[./perm_z]
type = PorousFlowPropertyAux
property = permeability
variable = perm_z
row = 2
column = 2
[../]
[./perm_x]
type = PorousFlowPropertyAux
property = permeability
variable = perm_x
row = 1
column = 1
[../]
[]
############################################################
[UserObjects]
[./fluid_mass_in_inc]
type = PorousFlowSumQuantity
[../]
[./fluid_mass_out_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_in_inc]
type = PorousFlowSumQuantity
[../]
[./heat_enthalpy_out_inc]
type = PorousFlowSumQuantity
[../]
[]
############################################################
[ICs]
[./init_pp]
type = FunctionIC
function = ${pp_ini_bc}
variable = porepressure
[../]
[./init_temp]
type = ConstantIC
value = ${T_ini_bc}
variable = temperature
[../]
[]
###############################################################
#[BCs]
# [./P_drained]
# type = DirichletBC
# variable = porepressure
# boundary = 'top bottom'
# value = ${pp_ini_bc}
# [../]
# [./T_cont]
# type = DirichletBC
# variable = temperature
# boundary = 'top bottom'
# value = ${T_ini_bc}
# [../]
#[]
############################################################
[DiracKernels]
[./summer_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 453.15
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./summer_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./summer_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_injection_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_injection_T]
type = EnthalpySink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_in_inc
line_base = '1 -9.99 -9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
pressure = porepressure
T_in = 399.15
fp = tabulated_water
p_or_t_vals = '-1e9 1e9'
fluxes = '-${inj_ext_flux} -${inj_ext_flux}'
[../]
[./winter_production_P]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = porepressure
SumQuantityUO = fluid_mass_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[./winter_production_T]
type = PorousFlowPolyLineSink
fluid_phase = 0
variable = temperature
SumQuantityUO = heat_enthalpy_out_inc
line_base = '1 9.99 9.99 -10.8'
line_length = ${well_length}
line_direction = '0 0 1'
use_mobility = false
use_enthalpy = true
p_or_t_vals = '-1e9 1e9'
fluxes = '${inj_ext_flux} ${inj_ext_flux}'
[../]
[]
############################################################
[Controls]
[summer_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_summer
enable_objects = 'DiracKernels::summer_injection_P DiracKernels::summer_injection_T DiracKernels::summer_production_P DiracKernels::summer_production_T'
[]
[winter_injection]
type = ConditionalFunctionEnableControl
conditional_function = inj_function_winter
enable_objects = 'DiracKernels::winter_injection_P DiracKernels::winter_injection_T DiracKernels::winter_production_P DiracKernels::winter_production_T'
[]
[]
[Functions]
[./inj_function_summer]
type = ParsedFunction
expression = '(t/24/3600/365-floor(t/24/3600/365))<=0.25'
[../]
[./inj_function_winter]
type = ParsedFunction
expression = '(t/24/3600/365-floor(t/24/3600/365))>0.5 & (t/24/3600/365-floor(t/24/3600/365)) <=0.75'
[../]
[]
############################################################
[FluidProperties]
[./true_water]
type = Water97FluidProperties
[../]
[./tabulated_water]
type = TabulatedFluidProperties
fp = true_water
allow_fp_and_tabulation = true
temperature_min = 275
temperature_max = 600
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_file = water97_tabulated.csv
[../]
[]
############################################################
[Materials]
[./internal_energy_aquifer]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 920.0
density = 2600.0
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
[./thermal_conductivity_aquifer]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.51 0 0 0 2.51 0 0 0 2.51'
block = 'inj_well ext_well caps layer_1 layer_2 layer_3 layer_4 layer_5 layer_6 layer_7 layer_8 layer_9 layer_10 layer_11 layer_12 layer_13 layer_14 layer_15 layer_16 layer_17 layer_18 layer_19 layer_20'
[../]
# inj_well/ext_well are separate blocks only in the new mesh (the old mesh had
# no dedicated well blocks, since its DiracKernel wells spanned the whole
# formation thickness). Borrowed verbatim from the newboundary input so every
# block in the (new) mesh has a defined porosity/permeability.
[./porosity_well]
type = PorousFlowPorosityConst
block = 'inj_well ext_well'
porosity = 0.299
[../]
[./permeability_well]
type = PorousFlowPermeabilityConst
block = 'inj_well ext_well'
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-09'
[../]
[./porosity_caps]
type = PorousFlowPorosityConst
block = caps
porosity = 0.01
[../]
[./permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-18 0 0 0 1E-18 0 0 0 1E-18'
[../]
[./porosity_layer_1]
type = PorousFlowPorosity
block = layer_1
porosity_zero = 0.143
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_1]
type = PorousFlowPermeabilityConst
block = layer_1
permeability = '.8130E-15 0 0 0 .8130E-15 0 0 0 .3160E-15'
[../]
[./porosity_layer_2]
type = PorousFlowPorosity
block = layer_2
porosity_zero = 0.158
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_2]
type = PorousFlowPermeabilityConst
block = layer_2
permeability = '.3460E-14 0 0 0 .3460E-14 0 0 0 .1350E-14'
[../]
[./porosity_layer_3]
type = PorousFlowPorosity
block = layer_3
porosity_zero = 0.176
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_3]
type = PorousFlowPermeabilityConst
block = layer_3
permeability = '.1120E-13 0 0 0 .1120E-13 0 0 0 .4390E-14'
[../]
[./porosity_layer_4]
type = PorousFlowPorosity
block = layer_4
porosity_zero = 0.188
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_4]
type = PorousFlowPermeabilityConst
block = layer_4
permeability = '.1920E-13 0 0 0 .1920E-13 0 0 0 .1400E-13'
[../]
[./porosity_layer_5]
type = PorousFlowPorosity
block = layer_5
porosity_zero = 0.166
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_5]
type = PorousFlowPermeabilityConst
block = layer_5
permeability = '.5500E-14 0 0 0 .5500E-14 0 0 0 .3900E-14'
[../]
[./porosity_layer_6]
type = PorousFlowPorosity
block = layer_6
porosity_zero = 0.185
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_6]
type = PorousFlowPermeabilityConst
block = layer_6
permeability = '.1770E-13 0 0 0 .1770E-13 0 0 0 .1190E-13'
[../]
[./porosity_layer_7]
type = PorousFlowPorosity
block = layer_7
porosity_zero = 0.24
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_7]
type = PorousFlowPermeabilityConst
block = layer_7
permeability = '.9110E-13 0 0 0 .9110E-13 0 0 0 .8870E-13'
[../]
[./porosity_layer_8]
type = PorousFlowPorosity
block = layer_8
porosity_zero = 0.26
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_8]
type = PorousFlowPermeabilityConst
block = layer_8
permeability = '.1390E-12 0 0 0 .1390E-12 0 0 0 .1290E-12'
[../]
[./porosity_layer_9]
type = PorousFlowPorosity
block = layer_9
porosity_zero = 0.285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_9]
type = PorousFlowPermeabilityConst
block = layer_9
permeability = '.2120E-12 0 0 0 .2120E-12 0 0 0 .2000E-12'
[../]
[./porosity_layer_10]
type = PorousFlowPorosity
block = layer_10
porosity_zero = 0.27
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_10]
type = PorousFlowPermeabilityConst
block = layer_10
permeability = '.1630E-12 0 0 0 .1630E-12 0 0 0 .1580E-12'
[../]
[./porosity_layer_11]
type = PorousFlowPorosity
block = layer_11
porosity_zero = 0.23
[../]
[./permeability_layer_11]
type = PorousFlowPermeabilityConst
block = layer_11
permeability = '.7700E-13 0 0 0 .7700E-13 0 0 0 .5990E-13'
[../]
[./porosity_layer_12]
type = PorousFlowPorosity
block = layer_12
porosity_zero = 0.306
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_12]
type = PorousFlowPermeabilityConst
block = layer_12
permeability = '.2880E-12 0 0 0 .2880E-12 0 0 0 .2880E-12'
[../]
[./porosity_layer_13]
type = PorousFlowPorosity
block = layer_13
porosity_zero = .285
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_13]
type = PorousFlowPermeabilityConst
block = layer_13
permeability = '.2130E-12 0 0 0 .2130E-12 0 0 0 .2040E-12'
[../]
[./porosity_layer_14]
type = PorousFlowPorosity
block = layer_14
porosity_zero = .155
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_14]
type = PorousFlowPermeabilityConst
block = layer_14
permeability = '.1080E-13 0 0 0 .1080E-13 0 0 0 .2190E-15'
[../]
[./porosity_layer_15]
type = PorousFlowPorosity
block = layer_15
porosity_zero = .299
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_15]
type = PorousFlowPermeabilityConst
block = layer_15
permeability = '.2660E-12 0 0 0 .2660E-12 0 0 0 .2330E-12'
[../]
[./porosity_layer_16]
type = PorousFlowPorosity
block = layer_16
porosity_zero = .334
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_16]
type = PorousFlowPermeabilityConst
block = layer_16
permeability = '.4230E-12 0 0 0 .4230E-12 0 0 0 .4110E-12'
[../]
[./porosity_layer_17]
type = PorousFlowPorosity
block = layer_17
porosity_zero = .32
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_17]
type = PorousFlowPermeabilityConst
block = layer_17
permeability = '.3560E-12 0 0 0 .3560E-12 0 0 0 .3300E-12'
[../]
[./porosity_layer_18]
type = PorousFlowPorosity
block = layer_18
porosity_zero = .294
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_18]
type = PorousFlowPermeabilityConst
block = layer_18
permeability = '.2420E-12 0 0 0 .2420E-12 0 0 0 .2360E-12'
[../]
[./porosity_layer_19]
type = PorousFlowPorosity
block = layer_19
porosity_zero = .279
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_19]
type = PorousFlowPermeabilityConst
block = layer_19
permeability = '.2050E-12 0 0 0 .2050E-12 0 0 0 .1310E-12'
[../]
[./porosity_layer_20]
type = PorousFlowPorosity
block = layer_20
porosity_zero = .071
fluid = true
biot_coefficient = 0.8
reference_porepressure = ${pp_ini_bc}
solid_bulk = 1.5E8
[../]
[./permeability_layer_20]
type = PorousFlowPermeabilityConst
block = layer_20
permeability = '.2700E-16 0 0 0 .2700E-16 0 0 0 .1000E-17 '
[../]
[]
############################################################
[Postprocessors]
[./is_summer]
type = FunctionValuePostprocessor
function = inj_function_summer
[../]
[./is_winter]
type = FunctionValuePostprocessor
function = inj_function_winter
[../]
[./inlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_in_inc
[../]
[./inlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_in_inc
[../]
[./outlet_mass_kg]
type = PorousFlowPlotQuantity
uo = fluid_mass_out_inc
[../]
[./outlet_enthalpy_J]
type = PorousFlowPlotQuantity
uo = heat_enthalpy_out_inc
[../]
[./production_temperature]
type = AverageNodalVariableValue
boundary = extraction
variable = temperature
[../]
[./production_porepressure]
type = AverageNodalVariableValue
boundary = extraction
variable = porepressure
[../]
[./injection_temperature]
type = AverageNodalVariableValue
boundary = injection
variable = temperature
[../]
[./injection_porepressure]
type = AverageNodalVariableValue
boundary = injection
variable = porepressure
[../]
[total]
type = MemoryUsage
mem_units = 'bytes'
execute_on = 'INITIAL TIMESTEP_END'
outputs = none # memory usage is not reproducible; excluded from regression CSV
[]
[]
############################################################
[Preconditioning]
active = 'ilu_may_use_less_mem'
#active = 'superlu'
[./ilu_may_use_less_mem]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels'
petsc_options_value = 'gmres asm ilu NONZERO 2'
[../]
[./original]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm lu NONZERO'
[../]
[./superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'gmres lu superlu_dist'
[../]
[./sub_pc_superlu]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_mat_solver_package'
petsc_options_value = 'gmres asm lu NONZERO superlu_dist'
[../]
[./preferred]
#per Andy
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[../]
[./basic]
type = SMP
full = true
petsc_options = '-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'
[../]
[]
##########################################################
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 31635000
# steady_state_detection = true
# steady_state_start_time = 864000 # 10 days
# steady_state_tolerance = 1e-6
dtmax = 216000 # 2 days
dtmin = 10
[./TimeStepper]
type = IterationAdaptiveDT
dt = 1000.0
optimal_iterations = 10
[../]
# controls for linear iterations
l_max_its = 100
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 50
nl_abs_tol = 1e-5
[]
###########################################################
[Outputs]
[./exodus]
type = Exodus
# file_base = sequence_out_1
# elemental_as_nodal = true
# overwrite = true
# sequence = true
time_step_interval = 20
[../]
[./Console]
type = Console
#perf_log = false
#output_linear = false
output_nonlinear = true
time_step_interval = 1
[../]
[./CSV]
type = CSV
time_step_interval = 1
[../]
[]
###########################################################