- 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
- boundaryThe list of boundary IDs from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundary IDs from the mesh where this object applies
- multipliersTuple of multiplying values. The flux values are multiplied by these.
C++ Type:std::vector<double>
Unit:(no unit assumed)
Controllable:No
Description:Tuple of multiplying values. The flux values are multiplied by these.
- pt_valsTuple of pressure values (for the fluid_phase specified). Must be monotonically increasing. For heat fluxes that don't involve fluids, these are temperature values
C++ Type:std::vector<double>
Unit:(no unit assumed)
Controllable:No
Description:Tuple of pressure values (for the fluid_phase specified). Must be monotonically increasing. For heat fluxes that don't involve fluids, these are temperature values
- 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
PorousFlowPiecewiseLinearSink
Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.
The basic sink is multiplied by a piecewise linear MOOSE Function of the pressure of a fluid phase or the temperature : Here the units of are kg.m.s (for fluids) or J.m.s (for heat).
If then the boundary condition will act as a sink, while if the boundary condition acts as a source. If applied to a fluid-component equation, the function has units kg.m.s. If applied to the heat equation, the function has units J.m.s. These units are potentially modified if the extra building blocks enumerated below are used.
In addition, the sink may be multiplied by any or all of the following quantities through the optional parameters list.
Fluid relative permeability
Fluid mobility (, where is the normal vector to the boundary)
Fluid mass fraction
Fluid internal energy
Thermal conductivity
See boundary conditions for many more details and discussion.
Input Parameters
- PT_shiftWhenever the sink is an explicit function of porepressure (such as a PiecewiseLinear function) the argument of the function is set to P - PT_shift instead of simply P. Similarly for temperature. PT_shift does not enter into any use_* calculations.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Whenever the sink is an explicit function of porepressure (such as a PiecewiseLinear function) the argument of the function is set to P - PT_shift instead of simply P. Similarly for temperature. PT_shift does not enter into any use_* calculations.
- displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
- fluid_phaseIf supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy. If not supplied, then this BC can only be a function of temperature
C++ Type:unsigned int
Controllable:No
Description:If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy. If not supplied, then this BC can only be a function of temperature
- flux_function1The flux. The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE. The functional form is useful for spatially or temporally varying sinks. Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)
Default:1
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:The flux. The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE. The functional form is useful for spatially or temporally varying sinks. Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)
- 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
- 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)
- use_enthalpyFalseIf true, then fluxes are multiplied by enthalpy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*
Default:False
C++ Type:bool
Controllable:No
Description:If true, then fluxes are multiplied by enthalpy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*
- use_internal_energyFalseIf true, then fluxes are multiplied by fluid internal energy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*
Default:False
C++ Type:bool
Controllable:No
Description:If true, then fluxes are multiplied by fluid internal energy. In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). This can be used in conjunction with other use_*
- use_mobilityFalseIf true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary. In this case bare_flux is measured in Pa.m^-1. This can be used in conjunction with other use_*
Default:False
C++ Type:bool
Controllable:No
Description:If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary. In this case bare_flux is measured in Pa.m^-1. This can be used in conjunction with other use_*
- use_relpermFalseIf true, then fluxes are multiplied by relative permeability. This can be used in conjunction with other use_*
Default:False
C++ Type:bool
Controllable:No
Description:If true, then fluxes are multiplied by relative permeability. This can be used in conjunction with other use_*
- use_thermal_conductivityFalseIf true, then fluxes are multiplied by thermal conductivity projected onto the normal direction. This can be used in conjunction with other use_*
Default:False
C++ Type:bool
Controllable:No
Description:If true, then fluxes are multiplied by thermal conductivity projected onto the normal direction. This can be used in conjunction with other use_*
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_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.
- diag_save_inThe name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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
- save_inThe name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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
- skip_execution_outside_variable_domainFalseWhether to skip execution of this boundary condition when the variable it applies to is not defined on the boundary. This can facilitate setups with moving variable domains and fixed boundaries. Note that the FEProblem boundary-restricted integrity checks will also need to be turned off if using this option
Default:False
C++ Type:bool
Controllable:No
Description:Whether to skip execution of this boundary condition when the variable it applies to is not defined on the boundary. This can facilitate setups with moving variable domains and fixed boundaries. Note that the FEProblem boundary-restricted integrity checks will also need to be turned off if using this option
- 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
- (modules/porous_flow/test/tests/jacobian/pls02.i)
- (modules/porous_flow/test/tests/newton_cooling/nc01.i)
- (modules/porous_flow/test/tests/numerical_diffusion/pffltvd.i)
- (modules/porous_flow/examples/tutorial/11_2D.i)
- (modules/porous_flow/test/tests/actions/basicthm_h.i)
- (modules/porous_flow/test/tests/jacobian/pls04.i)
- (modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_adaptivity.i)
- (modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)
- (modules/porous_flow/test/tests/sinks/s09.i)
- (modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_2D.i)
- (modules/porous_flow/test/tests/jacobian/pls03.i)
- (modules/porous_flow/test/tests/jacobian/pls01.i)
- (modules/porous_flow/test/tests/recover/theis.i)
- (modules/porous_flow/test/tests/actions/basicthm_hm.i)
- (modules/porous_flow/test/tests/sinks/s04.i)
- (modules/porous_flow/test/tests/actions/basicthm_thm.i)
- (modules/porous_flow/test/tests/actions/basicthm_th.i)
- (modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D_adaptivity.i)
- (modules/porous_flow/examples/restart/gas_injection_new_mesh.i)
- (modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D.i)
- (modules/porous_flow/test/tests/sinks/s09_fully_saturated.i)
- (modules/porous_flow/examples/restart/gas_injection.i)
- (modules/porous_flow/test/tests/sinks/PorousFlowPiecewiseLinearSink_BC_eg1.i)
- (modules/porous_flow/examples/co2_intercomparison/1Dradial/1Dradial.i)
- (modules/porous_flow/examples/tutorial/11.i)
- (modules/porous_flow/test/tests/newton_cooling/nc02.i)
- (modules/porous_flow/test/tests/recover/pffltvd.i)
- (modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_3D.i)
- (modules/porous_flow/test/tests/numerical_diffusion/no_action.i)
- (modules/porous_flow/test/tests/newton_cooling/nc04.i)
- (modules/porous_flow/test/tests/numerical_diffusion/pffltvd_action.i)
- (modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_2D_trimesh.i)
- (modules/porous_flow/examples/coal_mining/fine_with_fluid.i)
- (modules/porous_flow/test/tests/newton_cooling/nc06.i)
- (modules/porous_flow/test/tests/sinks/injection_production_eg.i)
- (modules/porous_flow/test/tests/numerical_diffusion/fully_saturated_action.i)
- (modules/porous_flow/test/tests/newton_cooling/nc08.i)
(modules/porous_flow/test/tests/jacobian/pls02.i)
# PorousFlowPiecewiseLinearSink with 2-phase, 2-components
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 3
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[ppwater]
[]
[ppgas]
[]
[massfrac_ph0_sp0]
[]
[massfrac_ph1_sp0]
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'ppwater ppgas massfrac_ph0_sp0 massfrac_ph1_sp0'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[ICs]
[ppwater]
type = RandomIC
variable = ppwater
min = -1
max = 0
[]
[ppgas]
type = RandomIC
variable = ppgas
min = 0
max = 1
[]
[massfrac_ph0_sp0]
type = RandomIC
variable = massfrac_ph0_sp0
min = 0
max = 1
[]
[massfrac_ph1_sp0]
type = RandomIC
variable = massfrac_ph1_sp0
min = 0
max = 1
[]
[]
[Kernels]
[dummy_ppwater]
type = TimeDerivative
variable = ppwater
[]
[dummy_ppgas]
type = TimeDerivative
variable = ppgas
[]
[dummy_m00]
type = TimeDerivative
variable = massfrac_ph0_sp0
[]
[dummy_m10]
type = TimeDerivative
variable = massfrac_ph1_sp0
[]
[]
[FluidProperties]
[simple_fluid0]
type = SimpleFluidProperties
bulk_modulus = 1.5
density0 = 1
thermal_expansion = 0
viscosity = 1
[]
[simple_fluid1]
type = SimpleFluidProperties
bulk_modulus = 0.5
density0 = 0.5
thermal_expansion = 0
viscosity = 1.4
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow2PhasePP
phase0_porepressure = ppwater
phase1_porepressure = ppgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
[]
[simple_fluid0]
type = PorousFlowSingleComponentFluid
fp = simple_fluid0
phase = 0
[]
[simple_fluid1]
type = PorousFlowSingleComponentFluid
fp = simple_fluid1
phase = 1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1 0 0 0 2 0 0 0 3'
[]
[relperm0]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[relperm1]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
[]
[]
[BCs]
[flux_w]
type = PorousFlowPiecewiseLinearSink
boundary = 'left'
pt_vals = '-1 -0.5 0'
multipliers = '1 2 4'
variable = ppwater
mass_fraction_component = 0
fluid_phase = 0
use_relperm = true
use_mobility = true
flux_function = 'x*y'
[]
[flux_g]
type = PorousFlowPiecewiseLinearSink
boundary = 'top'
pt_vals = '0 0.5 1'
multipliers = '1 -2 4'
mass_fraction_component = 0
variable = ppgas
fluid_phase = 1
use_relperm = true
use_mobility = true
flux_function = '-x*y'
[]
[flux_1]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '0 0.5 1'
multipliers = '1 3 4'
mass_fraction_component = 1
variable = massfrac_ph0_sp0
fluid_phase = 0
use_relperm = true
use_mobility = true
[]
[flux_2]
type = PorousFlowPiecewiseLinearSink
boundary = 'back top'
pt_vals = '0 0.5 1'
multipliers = '0 1 -3'
mass_fraction_component = 1
variable = massfrac_ph1_sp0
fluid_phase = 1
use_relperm = true
use_mobility = true
flux_function = '0.5*x*y'
[]
[]
[Preconditioning]
[check]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 2
[]
[Outputs]
file_base = pls02
[]
(modules/porous_flow/test/tests/newton_cooling/nc01.i)
# Newton cooling from a bar. 1-phase transient
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1000
ny = 1
xmin = 0
xmax = 100
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pressure'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Variables]
[pressure]
initial_condition = 2E6
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pressure
[]
[flux]
type = PorousFlowAdvectiveFlux
fluid_component = 0
gravity = '0 0 0'
variable = pressure
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1e6
density0 = 1000
thermal_expansion = 0
viscosity = 1e-3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
n = 2
phase = 0
[]
[]
[BCs]
[left]
type = DirichletBC
variable = pressure
boundary = left
value = 2E6
[]
[newton]
type = PorousFlowPiecewiseLinearSink
variable = pressure
boundary = right
pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 1
[]
[]
[VectorPostprocessors]
[porepressure]
type = LineValueSampler
variable = pressure
start_point = '0 0.5 0'
end_point = '100 0.5 0'
sort_by = x
num_points = 20
execute_on = timestep_end
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-12 1E-15 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E8
dt = 1E6
[]
[Outputs]
file_base = nc01
[along_line]
type = CSV
execute_vector_postprocessors_on = final
[]
[]
(modules/porous_flow/test/tests/numerical_diffusion/pffltvd.i)
# Using flux-limited TVD advection ala Kuzmin and Turek, employing PorousFlow Kernels and UserObjects, with superbee flux-limiter
[Mesh]
type = GeneratedMesh
dim = 1
nx = 100
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowFluxLimitedTVDAdvection
variable = tracer
advective_flux_calculator = advective_flux_calculator_0
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowFluxLimitedTVDAdvection
variable = porepressure
advective_flux_calculator = advective_flux_calculator_1
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[advective_flux_calculator_0]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 0
[]
[advective_flux_calculator_1]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 101
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-2
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/examples/tutorial/11_2D.i)
# Two-phase borehole injection problem in RZ coordinates
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
xmin = 1.0
xmax = 10
bias_x = 1.4
ny = 3
ymin = -6
ymax = 6
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 -2 0'
top_right = '10 2 0'
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x<1.0001'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
coord_type = RZ
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pwater pgas T disp_r'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[]
[]
[GlobalParams]
displacements = 'disp_r disp_z'
gravity = '0 0 0'
biot_coefficient = 1.0
PorousFlowDictator = dictator
[]
[Variables]
[pwater]
initial_condition = 20E6
[]
[pgas]
initial_condition = 20.1E6
[]
[T]
initial_condition = 330
scaling = 1E-5
[]
[disp_r]
scaling = 1E-5
[]
[]
[Kernels]
[mass_water_dot]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pwater
[]
[flux_water]
type = PorousFlowAdvectiveFlux
fluid_component = 0
use_displaced_mesh = false
variable = pwater
[]
[vol_strain_rate_water]
type = PorousFlowMassVolumetricExpansion
fluid_component = 0
variable = pwater
[]
[mass_co2_dot]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pgas
[]
[flux_co2]
type = PorousFlowAdvectiveFlux
fluid_component = 1
use_displaced_mesh = false
variable = pgas
[]
[vol_strain_rate_co2]
type = PorousFlowMassVolumetricExpansion
fluid_component = 1
variable = pgas
[]
[energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = T
[]
[advection]
type = PorousFlowHeatAdvection
use_displaced_mesh = false
variable = T
[]
[conduction]
type = PorousFlowHeatConduction
use_displaced_mesh = false
variable = T
[]
[vol_strain_rate_heat]
type = PorousFlowHeatVolumetricExpansion
variable = T
[]
[grad_stress_r]
type = StressDivergenceRZTensors
temperature = T
variable = disp_r
eigenstrain_names = thermal_contribution
use_displaced_mesh = false
component = 0
[]
[poro_r]
type = PorousFlowEffectiveStressCoupling
variable = disp_r
use_displaced_mesh = false
component = 0
[]
[]
[AuxVariables]
[disp_z]
[]
[effective_fluid_pressure]
family = MONOMIAL
order = CONSTANT
[]
[mass_frac_phase0_species0]
initial_condition = 1 # all water in phase=0
[]
[mass_frac_phase1_species0]
initial_condition = 0 # no water in phase=1
[]
[sgas]
family = MONOMIAL
order = CONSTANT
[]
[swater]
family = MONOMIAL
order = CONSTANT
[]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_tt]
family = MONOMIAL
order = CONSTANT
[]
[stress_zz]
family = MONOMIAL
order = CONSTANT
[]
[porosity]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[effective_fluid_pressure]
type = ParsedAux
coupled_variables = 'pwater pgas swater sgas'
expression = 'pwater * swater + pgas * sgas'
variable = effective_fluid_pressure
[]
[swater]
type = PorousFlowPropertyAux
variable = swater
property = saturation
phase = 0
execute_on = timestep_end
[]
[sgas]
type = PorousFlowPropertyAux
variable = sgas
property = saturation
phase = 1
execute_on = timestep_end
[]
[stress_rr_aux]
type = RankTwoAux
variable = stress_rr
rank_two_tensor = stress
index_i = 0
index_j = 0
[]
[stress_tt]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_tt
index_i = 2
index_j = 2
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 1
index_j = 1
[]
[porosity]
type = PorousFlowPropertyAux
variable = porosity
property = porosity
execute_on = timestep_end
[]
[]
[BCs]
[pinned_top_bottom_r]
type = DirichletBC
variable = disp_r
value = 0
boundary = 'top bottom'
[]
[cavity_pressure_r]
type = Pressure
boundary = injection_area
variable = disp_r
postprocessor = constrained_effective_fluid_pressure_at_wellbore
use_displaced_mesh = false
[]
[cold_co2]
type = DirichletBC
boundary = injection_area
variable = T
value = 290 # injection temperature
use_displaced_mesh = false
[]
[constant_co2_injection]
type = PorousFlowSink
boundary = injection_area
variable = pgas
fluid_phase = 1
flux_function = -1E-4
use_displaced_mesh = false
[]
[outer_water_removal]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pwater
fluid_phase = 0
pt_vals = '0 1E9'
multipliers = '0 1E8'
PT_shift = 20E6
use_mobility = true
use_relperm = true
use_displaced_mesh = false
[]
[outer_co2_removal]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pgas
fluid_phase = 1
pt_vals = '0 1E9'
multipliers = '0 1E8'
PT_shift = 20.1E6
use_mobility = true
use_relperm = true
use_displaced_mesh = false
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[tabulated_water]
type = TabulatedBicubicFluidProperties
fp = true_water
temperature_min = 275
pressure_max = 1E8
fluid_property_output_file = water97_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated_11.csv
[]
[true_co2]
type = CO2FluidProperties
[]
[tabulated_co2]
type = TabulatedBicubicFluidProperties
fp = true_co2
temperature_min = 275
pressure_max = 1E8
fluid_property_output_file = co2_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = co2_tabulated_11.csv
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = T
[]
[saturation_calculator]
type = PorousFlow2PhasePP
phase0_porepressure = pwater
phase1_porepressure = pgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
[]
[water]
type = PorousFlowSingleComponentFluid
fp = tabulated_water
phase = 0
[]
[co2]
type = PorousFlowSingleComponentFluid
fp = tabulated_co2
phase = 1
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 4
s_res = 0.1
sum_s_res = 0.2
phase = 0
[]
[relperm_co2]
type = PorousFlowRelativePermeabilityBC
nw_phase = true
lambda = 2
s_res = 0.1
sum_s_res = 0.2
phase = 1
[]
[porosity]
type = PorousFlowPorosity
fluid = true
mechanical = true
thermal = true
porosity_zero = 0.1
reference_temperature = 330
reference_porepressure = 20E6
thermal_expansion_coeff = 15E-6 # volumetric
solid_bulk = 8E9 # unimportant since biot = 1
[]
[permeability_aquifer]
type = PorousFlowPermeabilityKozenyCarman
block = aquifer
poroperm_function = kozeny_carman_phi0
phi0 = 0.1
n = 2
m = 2
k0 = 1E-12
[]
[permeability_caps]
type = PorousFlowPermeabilityKozenyCarman
block = caps
poroperm_function = kozeny_carman_phi0
phi0 = 0.1
n = 2
m = 2
k0 = 1E-15
k_anisotropy = '1 0 0 0 1 0 0 0 0.1'
[]
[rock_thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2 0 0 0 2 0 0 0 2'
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 1100
density = 2300
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5E9
poissons_ratio = 0.0
[]
[strain]
type = ComputeAxisymmetricRZSmallStrain
eigenstrain_names = 'thermal_contribution initial_stress'
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = T
thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
eigenstrain_name = thermal_contribution
stress_free_temperature = 330
[]
[initial_strain]
type = ComputeEigenstrainFromInitialStress
initial_stress = '20E6 0 0 0 20E6 0 0 0 20E6'
eigenstrain_name = initial_stress
[]
[stress]
type = ComputeLinearElasticStress
[]
[effective_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[volumetric_strain]
type = PorousFlowVolumetricStrain
[]
[]
[Postprocessors]
[effective_fluid_pressure_at_wellbore]
type = PointValue
variable = effective_fluid_pressure
point = '1 0 0'
execute_on = timestep_begin
use_displaced_mesh = false
[]
[constrained_effective_fluid_pressure_at_wellbore]
type = FunctionValuePostprocessor
function = constrain_effective_fluid_pressure
execute_on = timestep_begin
[]
[]
[Functions]
[constrain_effective_fluid_pressure]
type = ParsedFunction
symbol_names = effective_fluid_pressure_at_wellbore
symbol_values = effective_fluid_pressure_at_wellbore
expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E3
[TimeStepper]
type = IterationAdaptiveDT
dt = 1E3
growth_factor = 1.2
optimal_iterations = 10
[]
nl_abs_tol = 1E-7
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/actions/basicthm_h.i)
# PorousFlowBasicTHM action with coupling_type = HydroGenerator
# (no thermal or mechanical effects)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[]
[AuxVariables]
[temperature]
initial_condition = 293
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = Hydro
gravity = '0 0 0'
fp = simple_fluid
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-15
nl_rel_tol = 1E-14
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/jacobian/pls04.i)
# PorousFlowPiecewiseLinearSink with 2-phase, 3-components, with enthalpy, internal_energy, and thermal_conductivity
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 2
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[ppwater]
[]
[ppgas]
[]
[massfrac_ph0_sp0]
[]
[massfrac_ph0_sp1]
[]
[massfrac_ph1_sp0]
[]
[massfrac_ph1_sp1]
[]
[temp]
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp ppwater ppgas massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
number_fluid_phases = 2
number_fluid_components = 3
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[ICs]
[temp]
type = RandomIC
variable = temp
min = 1
max = 2
[]
[ppwater]
type = RandomIC
variable = ppwater
min = -1
max = 0
[]
[ppgas]
type = RandomIC
variable = ppgas
min = 0
max = 1
[]
[massfrac_ph0_sp0]
type = RandomIC
variable = massfrac_ph0_sp0
min = 0
max = 1
[]
[massfrac_ph0_sp1]
type = RandomIC
variable = massfrac_ph0_sp1
min = 0
max = 1
[]
[massfrac_ph1_sp0]
type = RandomIC
variable = massfrac_ph1_sp0
min = 0
max = 1
[]
[massfrac_ph1_sp1]
type = RandomIC
variable = massfrac_ph1_sp1
min = 0
max = 1
[]
[]
[Kernels]
[dummy_temp]
type = TimeDerivative
variable = temp
[]
[dummy_ppwater]
type = TimeDerivative
variable = ppwater
[]
[dummy_ppgas]
type = TimeDerivative
variable = ppgas
[]
[dummy_m00]
type = TimeDerivative
variable = massfrac_ph0_sp0
[]
[dummy_m01]
type = TimeDerivative
variable = massfrac_ph0_sp1
[]
[dummy_m10]
type = TimeDerivative
variable = massfrac_ph1_sp0
[]
[dummy_m11]
type = TimeDerivative
variable = massfrac_ph1_sp1
[]
[]
[FluidProperties]
[simple_fluid0]
type = SimpleFluidProperties
bulk_modulus = 1.5
density0 = 1
thermal_expansion = 0
viscosity = 1
cv = 1.1
[]
[simple_fluid1]
type = SimpleFluidProperties
bulk_modulus = 0.5
density0 = 0.5
thermal_expansion = 0
viscosity = 1.4
cv = 1.8
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[ppss]
type = PorousFlow2PhasePP
phase0_porepressure = ppwater
phase1_porepressure = ppgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
[]
[simple_fluid0]
type = PorousFlowSingleComponentFluid
fp = simple_fluid0
phase = 0
[]
[simple_fluid1]
type = PorousFlowSingleComponentFluid
fp = simple_fluid1
phase = 1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1 0 0 0 2 0 0 0 3'
[]
[relperm0]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[relperm1]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '0.1 0.2 0.3 0.2 0 0.1 0.3 0.1 0.1'
wet_thermal_conductivity = '10 2 31 2 40 1 31 1 10'
exponent = 0.5
[]
[]
[BCs]
[flux_w]
type = PorousFlowPiecewiseLinearSink
boundary = 'left'
pt_vals = '-1 -0.5 0'
multipliers = '1 2 4'
variable = ppwater
mass_fraction_component = 0
fluid_phase = 0
use_relperm = true
use_mobility = true
use_enthalpy = true
flux_function = 'x*y'
[]
[flux_g]
type = PorousFlowPiecewiseLinearSink
boundary = 'top'
pt_vals = '0 0.5 1'
multipliers = '1 -2 4'
mass_fraction_component = 0
variable = ppgas
fluid_phase = 1
use_relperm = true
use_mobility = true
use_internal_energy = true
flux_function = '-x*y'
[]
[flux_1]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '0 0.5 1'
multipliers = '1 3 4'
mass_fraction_component = 1
variable = massfrac_ph0_sp0
fluid_phase = 0
use_relperm = true
use_mobility = true
use_internal_energy = true
[]
[flux_2]
type = PorousFlowPiecewiseLinearSink
boundary = 'back top'
pt_vals = '0 0.5 1'
multipliers = '0 1 -3'
mass_fraction_component = 1
variable = massfrac_ph1_sp0
fluid_phase = 1
use_relperm = true
use_mobility = true
use_enthalpy = true
flux_function = '0.5*x*y'
[]
[flux_3]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '0 0.5 1'
multipliers = '1 3 4'
mass_fraction_component = 2
variable = ppwater
fluid_phase = 0
use_relperm = true
use_enthalpy = true
use_mobility = true
[]
[flux_4]
type = PorousFlowPiecewiseLinearSink
boundary = 'back top'
pt_vals = '0 0.5 1'
multipliers = '0 1 -3'
mass_fraction_component = 2
variable = massfrac_ph1_sp0
fluid_phase = 1
use_relperm = true
use_mobility = true
flux_function = '-0.5*x*y'
use_enthalpy = true
use_thermal_conductivity = true
[]
[]
[Preconditioning]
[check]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
file_base = pls04
[]
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_adaptivity.i)
# Pressure pulse in 1D with 1 phase - transient simulation with a constant
# PorousFlowPorosity and mesh adaptivity with an indicator
[Mesh]
type = GeneratedMesh
dim = 1
nx = 10
xmin = 0
xmax = 100
[]
[Adaptivity]
marker = marker
[Markers]
[marker]
type = ErrorFractionMarker
indicator = front
refine = 0.5
coarsen = 0.2
[]
[]
[Indicators]
[front]
type = GradientJumpIndicator
variable = pp
[]
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 2E6
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[flux]
type = PorousFlowAdvectiveFlux
variable = pp
gravity = '0 0 0'
fluid_component = 0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1e-7
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
density0 = 1000
thermal_expansion = 0
viscosity = 1e-3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0
phase = 0
[]
[]
[BCs]
[left]
type = DirichletBC
boundary = left
preset = false
value = 3E6
variable = pp
[]
[right]
type = PorousFlowPiecewiseLinearSink
variable = pp
boundary = right
fluid_phase = 0
pt_vals = '0 1E9'
multipliers = '0 1E9'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E-6
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1e3
end_time = 5e3
[]
[Postprocessors]
[p000]
type = PointValue
variable = pp
point = '0 0 0'
execute_on = 'initial timestep_end'
[]
[p010]
type = PointValue
variable = pp
point = '10 0 0'
execute_on = 'initial timestep_end'
[]
[p020]
type = PointValue
variable = pp
point = '20 0 0'
execute_on = 'initial timestep_end'
[]
[p030]
type = PointValue
variable = pp
point = '30 0 0'
execute_on = 'initial timestep_end'
[]
[p040]
type = PointValue
variable = pp
point = '40 0 0'
execute_on = 'initial timestep_end'
[]
[p050]
type = PointValue
variable = pp
point = '50 0 0'
execute_on = 'initial timestep_end'
[]
[p060]
type = PointValue
variable = pp
point = '60 0 0'
execute_on = 'initial timestep_end'
[]
[p070]
type = PointValue
variable = pp
point = '70 0 0'
execute_on = 'initial timestep_end'
[]
[p080]
type = PointValue
variable = pp
point = '80 0 0'
execute_on = 'initial timestep_end'
[]
[p090]
type = PointValue
variable = pp
point = '90 0 0'
execute_on = 'initial timestep_end'
[]
[p100]
type = PointValue
variable = pp
point = '100 0 0'
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
print_linear_residuals = false
csv = true
[]
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)
# Strata deformation and fluid flow aaround a coal mine - 3D model
#
# A "half model" is used. The mine is 400m deep and
# just the roof is studied (-400<=z<=0). The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long. The outer boundaries
# are 1km from the excavation boundaries.
#
# The excavation takes 0.5 years.
#
# The boundary conditions for this simulation are:
# - disp_x = 0 at x=0 and x=1150
# - disp_y = 0 at y=-1000 and y=1000
# - disp_z = 0 at z=-400, but there is a time-dependent
# Young modulus that simulates excavation
# - wc_x = 0 at y=-1000 and y=1000
# - wc_y = 0 at x=0 and x=1150
# - no flow at x=0, z=-400 and z=0
# - fixed porepressure at y=-1000, y=1000 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# A single-phase unsaturated fluid is used.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa, and time units are measured in years.
#
# The initial porepressure is hydrostatic with P=0 at z=0, so
# Porepressure ~ - 0.01*z MPa, where the fluid has density 1E3 kg/m^3 and
# gravity = = 10 m.s^-2 = 1E-5 MPa m^2/kg.
# To be more accurate, i use
# Porepressure = -bulk * log(1 + g*rho0*z/bulk)
# where bulk=2E3 MPa and rho0=1Ee kg/m^3.
# The initial stress is consistent with the weight force from undrained
# density 2500 kg/m^3, and fluid porepressure, and a Biot coefficient of 0.7, ie,
# stress_zz^effective = 0.025*z + 0.7 * initial_porepressure
# The maximum and minimum principal horizontal effective stresses are
# assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 2 MPa
# MC friction angle = 35 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
# Fluid density at zero porepressure = 1E3 kg/m^3
# Fluid bulk modulus = 2E3 MPa
# Fluid viscosity = 1.1E-3 Pa.s = 1.1E-9 MPa.s = 3.5E-17 MPa.year
#
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
PorousFlowDictator = dictator
biot_coefficient = 0.7
[]
[Mesh]
[file]
type = FileMeshGenerator
file = mesh/coarse.e
[]
[xmin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmin
normal = '-1 0 0'
input = file
[]
[xmax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmax
normal = '1 0 0'
input = xmin
[]
[ymin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymin
normal = '0 -1 0'
input = xmax
[]
[ymax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymax
normal = '0 1 0'
input = ymin
[]
[zmax]
type = SideSetsAroundSubdomainGenerator
block = 16
new_boundary = zmax
normal = '0 0 1'
input = ymax
[]
[zmin]
type = SideSetsAroundSubdomainGenerator
block = 2
new_boundary = zmin
normal = '0 0 -1'
input = zmax
[]
[excav]
type = SubdomainBoundingBoxGenerator
input = zmin
block_id = 1
bottom_left = '0 0 -400'
top_right = '150 1000 -397'
[]
[roof]
type = SideSetsBetweenSubdomainsGenerator
primary_block = 3
paired_block = 1
input = excav
new_boundary = roof
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[wc_x]
[]
[wc_y]
[]
[porepressure]
scaling = 1E-5
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = ini_pp
[]
[]
[Kernels]
[cx_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_x
component = 0
[]
[cy_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_y
component = 1
[]
[cz_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_z
component = 2
[]
[x_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_x
displacements = 'wc_x wc_y wc_z'
component = 0
base_name = couple
[]
[y_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_y
displacements = 'wc_x wc_y wc_z'
component = 1
base_name = couple
[]
[x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[]
[y_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_y
component = 1
[]
[gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
component = 2
variable = disp_z
[]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[]
[flux]
type = PorousFlowAdvectiveFlux
use_displaced_mesh = false
variable = porepressure
gravity = '0 0 -10E-6'
fluid_component = 0
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
variable = porepressure
fluid_component = 0
[]
[]
[AuxVariables]
[saturation]
order = CONSTANT
family = MONOMIAL
[]
[darcy_x]
order = CONSTANT
family = MONOMIAL
[]
[darcy_y]
order = CONSTANT
family = MONOMIAL
[]
[darcy_z]
order = CONSTANT
family = MONOMIAL
[]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[wc_z]
[]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[stress_yx]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_zx]
order = CONSTANT
family = MONOMIAL
[]
[stress_zy]
order = CONSTANT
family = MONOMIAL
[]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_xx]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_xy]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_xz]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_yx]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_yz]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_zx]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_zy]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_zz]
order = CONSTANT
family = MONOMIAL
[]
[perm_xx]
order = CONSTANT
family = MONOMIAL
[]
[perm_yy]
order = CONSTANT
family = MONOMIAL
[]
[perm_zz]
order = CONSTANT
family = MONOMIAL
[]
[mc_shear]
order = CONSTANT
family = MONOMIAL
[]
[mc_tensile]
order = CONSTANT
family = MONOMIAL
[]
[wp_shear]
order = CONSTANT
family = MONOMIAL
[]
[wp_tensile]
order = CONSTANT
family = MONOMIAL
[]
[wp_shear_f]
order = CONSTANT
family = MONOMIAL
[]
[wp_tensile_f]
order = CONSTANT
family = MONOMIAL
[]
[mc_shear_f]
order = CONSTANT
family = MONOMIAL
[]
[mc_tensile_f]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[saturation_water]
type = PorousFlowPropertyAux
variable = saturation
property = saturation
phase = 0
execute_on = timestep_end
[]
[darcy_x]
type = PorousFlowDarcyVelocityComponent
variable = darcy_x
gravity = '0 0 -10E-6'
component = x
[]
[darcy_y]
type = PorousFlowDarcyVelocityComponent
variable = darcy_y
gravity = '0 0 -10E-6'
component = y
[]
[darcy_z]
type = PorousFlowDarcyVelocityComponent
variable = darcy_z
gravity = '0 0 -10E-6'
component = z
[]
[porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
execute_on = timestep_end
[]
[stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[]
[stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[]
[stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[]
[stress_yx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[]
[stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[]
[stress_zx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[]
[stress_zy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[]
[total_strain_xx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[]
[total_strain_xy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[]
[total_strain_xz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[]
[total_strain_yx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[]
[total_strain_yy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[]
[total_strain_yz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[]
[total_strain_zx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[]
[total_strain_zy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[]
[total_strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[]
[perm_xx]
type = PorousFlowPropertyAux
property = permeability
variable = perm_xx
row = 0
column = 0
execute_on = timestep_end
[]
[perm_yy]
type = PorousFlowPropertyAux
property = permeability
variable = perm_yy
row = 1
column = 1
execute_on = timestep_end
[]
[perm_zz]
type = PorousFlowPropertyAux
property = permeability
variable = perm_zz
row = 2
column = 2
execute_on = timestep_end
[]
[mc_shear]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_internal_parameter
variable = mc_shear
execute_on = timestep_end
[]
[mc_tensile]
type = MaterialStdVectorAux
index = 1
property = mc_plastic_internal_parameter
variable = mc_tensile
execute_on = timestep_end
[]
[wp_shear]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_internal_parameter
variable = wp_shear
execute_on = timestep_end
[]
[wp_tensile]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_internal_parameter
variable = wp_tensile
execute_on = timestep_end
[]
[mc_shear_f]
type = MaterialStdVectorAux
index = 6
property = mc_plastic_yield_function
variable = mc_shear_f
execute_on = timestep_end
[]
[mc_tensile_f]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_yield_function
variable = mc_tensile_f
execute_on = timestep_end
[]
[wp_shear_f]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_yield_function
variable = wp_shear_f
execute_on = timestep_end
[]
[wp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_yield_function
variable = wp_tensile_f
execute_on = timestep_end
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_x
boundary = 'xmin xmax'
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_y
boundary = 'ymin ymax'
value = 0.0
[]
[no_z]
type = DirichletBC
variable = disp_z
boundary = zmin
value = 0.0
[]
[no_wc_x]
type = DirichletBC
variable = wc_x
boundary = 'ymin ymax'
value = 0.0
[]
[no_wc_y]
type = DirichletBC
variable = wc_y
boundary = 'xmin xmax'
value = 0.0
[]
[fix_porepressure]
type = FunctionDirichletBC
variable = porepressure
boundary = 'ymin ymax xmax'
function = ini_pp
[]
[roof_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
pt_vals = '-1E3 1E3'
multipliers = '-1 1'
fluid_phase = 0
flux_function = roof_conductance
boundary = roof
[]
[roof_bcs]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = roof
[]
[]
[Functions]
[ini_pp]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0'
symbol_values = '2E3 0.0 1E-5 1E3'
expression = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
[]
[ini_xx]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
[]
[ini_zz]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
[]
[excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '0.5 0 1000.0 1E-9 1 60'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
[]
[density_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '0.5 0 1000.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[]
[roof_conductance]
type = ParsedFunction
symbol_names = 'end_t ymin ymax maxval minval'
symbol_values = '0.5 0 1000.0 1E7 0'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),maxval,minval)'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1 # MPa^-1
[]
[mc_coh_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.99 # MPa
value_residual = 2.01 # MPa
rate = 1.0
[]
[mc_fric]
type = TensorMechanicsHardeningConstant
value = 0.61 # 35deg
[]
[mc_dil]
type = TensorMechanicsHardeningConstant
value = 0.15 # 8deg
[]
[mc_tensile_str_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[]
[mc_compressive_str]
type = TensorMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[]
[wp_coh_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[]
[wp_tan_fric]
type = TensorMechanicsHardeningConstant
value = 0.26 # 15deg
[]
[wp_tan_dil]
type = TensorMechanicsHardeningConstant
value = 0.18 # 10deg
[]
[wp_tensile_str_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[]
[wp_compressive_str_soften]
type = TensorMechanicsHardeningCubic
value_0 = 100
value_residual = 1
internal_limit = 1.0
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E3
density0 = 1000
thermal_expansion = 0
viscosity = 3.5E-17
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity_bulk]
type = PorousFlowPorosity
fluid = true
mechanical = true
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
ensure_positive = true
porosity_zero = 0.02
solid_bulk = 5.3333E3
[]
[porosity_excav]
type = PorousFlowPorosityConst
block = 1
porosity = 1.0
[]
[permeability_bulk]
type = PorousFlowPermeabilityKozenyCarman
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
poroperm_function = kozeny_carman_phi0
k0 = 1E-15
phi0 = 0.02
n = 2
m = 2
[]
[permeability_excav]
type = PorousFlowPermeabilityConst
block = 1
permeability = '0 0 0 0 0 0 0 0 0'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 4
s_res = 0.4
sum_s_res = 0.4
phase = 0
[]
[elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
[]
[elasticity_tensor_1]
type = ComputeLayeredCosseratElasticityTensor
block = 1
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
elasticity_tensor_prefactor = excav_sideways
[]
[strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[]
[ini_stress]
type = ComputeEigenstrainFromInitialStress
eigenstrain_name = ini_stress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
[]
[stress_0]
type = ComputeMultipleInelasticCosseratStress
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
inelastic_models = 'mc wp'
cycle_models = true
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[]
[stress_1]
type = ComputeMultipleInelasticCosseratStress
block = 1
inelastic_models = ''
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[]
[mc]
type = CappedMohrCoulombCosseratStressUpdate
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = mc
tensile_strength = mc_tensile_str_strong_harden
compressive_strength = mc_compressive_str
cohesion = mc_coh_strong_harden
friction_angle = mc_fric
dilation_angle = mc_dil
max_NR_iterations = 100000
smoothing_tol = 0.1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-9 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[]
[wp]
type = CappedWeakPlaneCosseratStressUpdate
warn_about_precision_loss = false
base_name = wp
cohesion = wp_coh_harden
tan_friction_angle = wp_tan_fric
tan_dilation_angle = wp_tan_dil
tensile_strength = wp_tensile_str_harden
compressive_strength = wp_compressive_str_soften
max_NR_iterations = 10000
tip_smoother = 0.05
smoothing_tol = 0.05 # MPa # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0E-3
[]
[undrained_density_0]
type = GenericConstantMaterial
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
prop_names = density
prop_values = 2500
[]
[undrained_density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Postprocessors]
[min_roof_disp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = disp_z
[]
[min_roof_pp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = porepressure
[]
[min_surface_disp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = disp_z
[]
[min_surface_pp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = porepressure
[]
[max_perm_zz]
type = ElementExtremeValue
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
variable = perm_zz
[]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
# best overall
# petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
# petsc_options_value = ' lu mumps'
# best if you do not have mumps:
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu superlu_dist'
# best if you do not have mumps or superlu_dist:
#petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' asm 2 lu gmres 200'
# very basic:
#petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' bjacobi gmres 200'
line_search = bt
nl_abs_tol = 1e-3
nl_rel_tol = 1e-5
l_max_its = 200
nl_max_its = 30
start_time = 0.0
dt = 0.014706
end_time = 0.014706 #0.5
[]
[Outputs]
time_step_interval = 1
print_linear_residuals = true
exodus = true
csv = true
console = true
[]
(modules/porous_flow/test/tests/sinks/s09.i)
# Apply a piecewise-linear sink flux to the right-hand side and watch fluid flow to it
#
# This test has a single phase with two components. The test initialises with
# the porous material fully filled with component=1. The left-hand side is fixed
# at porepressure=1 and mass-fraction of the zeroth component being unity.
# The right-hand side has a very strong piecewise-linear flux that keeps the
# porepressure~0 at that side. Fluid mass is extracted by this flux in proportion
# to the fluid component mass fraction.
#
# Therefore, the zeroth fluid component will flow from left to right (down the
# pressure gradient).
#
# The important DE is
# porosity * dc/dt = (perm / visc) * grad(P) * grad(c)
# which is true for c = mass-fraction, and very large bulk modulus of the fluid.
# For grad(P) constant in time and space (as in this example) this is just the
# advection equation for c, with velocity = perm / visc / porosity. The parameters
# are chosen to velocity = 1 m/s.
# In the numerical world, and especially with full upwinding, the advection equation
# suffers from diffusion. In this example, the diffusion is obvious when plotting
# the mass-fraction along the line, but the average velocity of the front is still
# correct at 1 m/s.
[Mesh]
type = GeneratedMesh
dim = 1
nx = 100
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp frac'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[Variables]
[pp]
[]
[frac]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
function = 1-x
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = frac
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pp
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
gravity = '0 0 0'
variable = frac
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
gravity = '0 0 0'
variable = pp
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1e10 # need large in order for constant-velocity advection
density0 = 1 # almost irrelevant, except that the ability of the right BC to keep P fixed at zero is related to density_P0
thermal_expansion = 0
viscosity = 11
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = frac
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.1 0 0 0 1.1 0 0 0 1.1'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2 # irrelevant in this fully-saturated situation
phase = 0
[]
[]
[BCs]
[lhs_fixed_a]
type = DirichletBC
boundary = 'left'
variable = frac
value = 1
[]
[lhs_fixed_b]
type = DirichletBC
boundary = 'left'
variable = pp
value = 1
[]
[flux0]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '-100 100'
multipliers = '-1 1'
variable = frac # the zeroth comonent
mass_fraction_component = 0
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 1E4
[]
[flux1]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '-100 100'
multipliers = '-1 1'
variable = pp # comonent 1
mass_fraction_component = 1
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 1E4
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1E-2
end_time = 1
nl_rel_tol = 1E-12
nl_abs_tol = 1E-12
[]
[VectorPostprocessors]
[mf]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 100
sort_by = x
variable = frac
[]
[]
[Outputs]
file_base = s09
[console]
type = Console
execute_on = 'nonlinear linear'
[]
[csv]
type = CSV
sync_times = '0.1 0.5 1'
sync_only = true
[]
time_step_interval = 10
[]
(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_2D.i)
# Using flux-limited TVD advection ala Kuzmin and Turek, employing PorousFlow Kernels and UserObjects, with superbee flux-limiter
# 3D version
[Mesh]
type = GeneratedMesh
dim = 2
nx = 10
xmin = 0
xmax = 1
ny = 4
ymin = 0
ymax = 0.5
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowFluxLimitedTVDAdvection
variable = tracer
advective_flux_calculator = advective_flux_calculator_0
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowFluxLimitedTVDAdvection
variable = porepressure
advective_flux_calculator = advective_flux_calculator_1
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[advective_flux_calculator_0]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 0
[]
[advective_flux_calculator_1]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0.5 0'
num_points = 11
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-2
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/jacobian/pls03.i)
# PorousFlowPiecewiseLinearSink with 2-phase, 3-components
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 2
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[ppwater]
[]
[ppgas]
[]
[massfrac_ph0_sp0]
[]
[massfrac_ph0_sp1]
[]
[massfrac_ph1_sp0]
[]
[massfrac_ph1_sp1]
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'ppwater ppgas massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
number_fluid_phases = 2
number_fluid_components = 3
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[ICs]
[ppwater]
type = RandomIC
variable = ppwater
min = -1
max = 0
[]
[ppgas]
type = RandomIC
variable = ppgas
min = 0
max = 1
[]
[massfrac_ph0_sp0]
type = RandomIC
variable = massfrac_ph0_sp0
min = 0
max = 1
[]
[massfrac_ph0_sp1]
type = RandomIC
variable = massfrac_ph0_sp1
min = 0
max = 1
[]
[massfrac_ph1_sp0]
type = RandomIC
variable = massfrac_ph1_sp0
min = 0
max = 1
[]
[massfrac_ph1_sp1]
type = RandomIC
variable = massfrac_ph1_sp1
min = 0
max = 1
[]
[]
[Kernels]
[dummy_ppwater]
type = TimeDerivative
variable = ppwater
[]
[dummy_ppgas]
type = TimeDerivative
variable = ppgas
[]
[dummy_m00]
type = TimeDerivative
variable = massfrac_ph0_sp0
[]
[dummy_m01]
type = TimeDerivative
variable = massfrac_ph0_sp1
[]
[dummy_m10]
type = TimeDerivative
variable = massfrac_ph1_sp0
[]
[dummy_m11]
type = TimeDerivative
variable = massfrac_ph1_sp1
[]
[]
[FluidProperties]
[simple_fluid0]
type = SimpleFluidProperties
bulk_modulus = 1.5
density0 = 1
thermal_expansion = 0
viscosity = 1
[]
[simple_fluid1]
type = SimpleFluidProperties
bulk_modulus = 0.5
density0 = 0.5
thermal_expansion = 0
viscosity = 1.4
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow2PhasePP
phase0_porepressure = ppwater
phase1_porepressure = ppgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph0_sp1 massfrac_ph1_sp0 massfrac_ph1_sp1'
[]
[simple_fluid0]
type = PorousFlowSingleComponentFluid
fp = simple_fluid0
phase = 0
[]
[simple_fluid1]
type = PorousFlowSingleComponentFluid
fp = simple_fluid1
phase = 1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1 0 0 0 2 0 0 0 3'
[]
[relperm0]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[relperm1]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
[]
[]
[BCs]
[flux_w]
type = PorousFlowPiecewiseLinearSink
boundary = 'left'
pt_vals = '-1 -0.5 0'
multipliers = '1 2 4'
variable = ppwater
mass_fraction_component = 0
fluid_phase = 0
use_relperm = true
use_mobility = true
flux_function = 'x*y'
[]
[flux_g]
type = PorousFlowPiecewiseLinearSink
boundary = 'top'
pt_vals = '0 0.5 1'
multipliers = '1 -2 4'
mass_fraction_component = 0
variable = ppgas
fluid_phase = 1
use_relperm = true
use_mobility = true
flux_function = '-x*y'
[]
[flux_1]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '0 0.5 1'
multipliers = '1 3 4'
mass_fraction_component = 1
variable = massfrac_ph0_sp0
fluid_phase = 0
use_relperm = true
use_mobility = true
[]
[flux_2]
type = PorousFlowPiecewiseLinearSink
boundary = 'back top'
pt_vals = '0 0.5 1'
multipliers = '0 1 -3'
mass_fraction_component = 1
variable = massfrac_ph1_sp0
fluid_phase = 1
use_relperm = true
use_mobility = true
flux_function = '0.5*x*y'
[]
[flux_3]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '0 0.5 1'
multipliers = '1 3 4'
mass_fraction_component = 2
variable = ppwater
fluid_phase = 0
use_relperm = true
use_mobility = true
[]
[flux_4]
type = PorousFlowPiecewiseLinearSink
boundary = 'back top'
pt_vals = '0 0.5 1'
multipliers = '0 1 -3'
mass_fraction_component = 2
variable = massfrac_ph1_sp0
fluid_phase = 1
use_relperm = true
use_mobility = true
flux_function = '-0.5*x*y'
[]
[]
[Preconditioning]
[check]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 2
[]
[Outputs]
file_base = pls03
[]
(modules/porous_flow/test/tests/jacobian/pls01.i)
# PorousFlowPiecewiseLinearSink with 1-phase, 1-component
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[Variables]
[pp]
[]
[]
[ICs]
[pp]
type = RandomIC
variable = pp
max = 0
min = -1
[]
[]
[Kernels]
[dummy]
type = TimeDerivative
variable = pp
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1
density0 = 1
thermal_expansion = 0
viscosity = 1.1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.1 0 0 0 2.2 0 0 0 3.3'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[BCs]
[flux]
type = PorousFlowPiecewiseLinearSink
boundary = 'left'
pt_vals = '-1 -0.5 0'
multipliers = '1 2 4'
variable = pp
fluid_phase = 0
use_relperm = true
use_mobility = true
flux_function = 'x*y'
[]
[]
[Preconditioning]
[check]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 2
[]
[Outputs]
file_base = pls01
[]
(modules/porous_flow/test/tests/recover/theis.i)
# Tests that PorousFlow can successfully recover using a checkpoint file.
# This test contains stateful material properties, adaptivity and integrated
# boundary conditions with nodal-sized materials.
#
# This test file is run three times:
# 1) The full input file is run to completion
# 2) The input file is run for half the time and checkpointing is included
# 3) The input file is run in recovery using the checkpoint data
#
# The final output of test 3 is compared to the final output of test 1 to verify
# that recovery was successful.
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
nx = 20
xmax = 100
bias_x = 1.05
[]
coord_type = RZ
rz_coord_axis = Y
# To get consistent ordering of results with distributed meshes
allow_renumbering = false
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Adaptivity]
marker = marker
max_h_level = 4
[Indicators]
[front]
type = GradientJumpIndicator
variable = zi
[]
[]
[Markers]
[marker]
type = ErrorFractionMarker
indicator = front
refine = 0.8
coarsen = 0.2
[]
[]
[]
[AuxVariables]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[x1]
order = CONSTANT
family = MONOMIAL
[]
[y0]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = timestep_end
[]
[x1]
type = PorousFlowPropertyAux
variable = x1
property = mass_fraction
phase = 0
fluid_component = 1
execute_on = timestep_end
[]
[y0]
type = PorousFlowPropertyAux
variable = y0
property = mass_fraction
phase = 1
fluid_component = 0
execute_on = timestep_end
[]
[]
[Variables]
[pgas]
initial_condition = 20e6
[]
[zi]
initial_condition = 0
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pgas
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pgas
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = zi
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = zi
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
pc = 0
[]
[fs]
type = PorousFlowWaterNCG
water_fp = water
gas_fp = co2
capillary_pressure = pc
[]
[]
[FluidProperties]
[co2]
type = CO2FluidProperties
[]
[water]
type = Water97FluidProperties
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = 20
[]
[waterncg]
type = PorousFlowFluidState
gas_porepressure = pgas
z = zi
temperature_unit = Celsius
capillary_pressure = pc
fluid_state = fs
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.1
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 1
[]
[]
[BCs]
[aquifer]
type = PorousFlowPiecewiseLinearSink
variable = pgas
boundary = right
pt_vals = '0 1e8'
multipliers = '0 1e8'
flux_function = 1e-6
PT_shift = 20e6
[]
[]
[DiracKernels]
[source]
type = PorousFlowSquarePulsePointSource
point = '0 0 0'
mass_flux = 2
variable = zi
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 2e2
dt = 50
[]
[VectorPostprocessors]
[line]
type = NodalValueSampler
sort_by = x
variable = 'pgas zi'
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
csv = true
[]
(modules/porous_flow/test/tests/actions/basicthm_hm.i)
# PorousFlowBasicTHM action with coupling_type = HydroMechanical
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
displacements = 'disp_x disp_y'
biot_coefficient = 1.0
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[disp_x]
scaling = 1e-10
[]
[disp_y]
scaling = 1e-10
[]
[]
[AuxVariables]
[temperature]
initial_condition = 293
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = HydroMechanical
gravity = '0 0 0'
fp = simple_fluid
use_displaced_mesh = false
add_stress_aux = false
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[top_bottom]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'top bottom'
[]
[right]
type = DirichletBC
variable = disp_x
value = 0
boundary = right
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5e9
poissons_ratio = 0.0
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-14
nl_rel_tol = 1e-14
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/sinks/s04.i)
# apply a piecewise-linear sink flux and observe the correct behavior
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 1
zmin = 0
zmax = 2
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[Variables]
[pp]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
function = y+1
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1.3
density0 = 1.1
thermal_expansion = 0
viscosity = 1.1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-5 0 0 0 1E-5 0 0 0 1E-5'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[]
[AuxVariables]
[flux_out]
[]
[xval]
[]
[yval]
[]
[pt_shift]
initial_condition = 0.3
[]
[]
[ICs]
[xval]
type = FunctionIC
variable = xval
function = x
[]
[yval]
type = FunctionIC
variable = yval
function = y
[]
[]
[Postprocessors]
[p00]
type = PointValue
point = '0 0 0'
variable = pp
execute_on = 'initial timestep_end'
[]
[p10]
type = PointValue
point = '1 0 0'
variable = pp
execute_on = 'initial timestep_end'
[]
[m10]
type = ParsedPostprocessor
expression = 'vol*por*dens0*exp(p10/bulk)'
constant_names = 'vol por dens0 bulk'
constant_expressions = '0.25 0.1 1.1 1.3'
pp_names = p10
execute_on = 'initial timestep_end'
[]
[dm10]
type = ChangeOverTimePostprocessor
postprocessor = m10
outputs = none
[]
[m10_prev]
type = ParsedPostprocessor
expression = 'm10 - dm10'
pp_names = 'm10 dm10'
outputs = 'console'
[]
[m10_rate]
type = ParsedPostprocessor
expression = 'fcn*if(p10>0.8,1,if(p10<0.3,0.5,0.2+p10))'
constant_names = 'fcn'
constant_expressions = '8'
pp_names = 'p10'
[]
[m10_expect]
type = ParsedPostprocessor
expression = 'm10_prev-m10_rate*area*dt'
constant_names = 'area dt'
constant_expressions = '0.5 1E-3'
pp_names = 'm10_prev m10_rate'
[]
[p01]
type = PointValue
point = '0 1 0'
variable = pp
execute_on = 'initial timestep_end'
[]
[p11]
type = PointValue
point = '1 1 0'
variable = pp
execute_on = 'initial timestep_end'
[]
[m11]
type = ParsedPostprocessor
expression = 'vol*por*dens0*exp(p11/bulk)'
constant_names = 'vol por dens0 bulk'
constant_expressions = '0.25 0.1 1.1 1.3'
pp_names = 'p11'
execute_on = 'initial timestep_end'
[]
[dm11]
type = ChangeOverTimePostprocessor
postprocessor = m11
outputs = none
[]
[m11_prev]
type = ParsedPostprocessor
expression = 'm11 - dm11'
pp_names = 'm11 dm11'
outputs = 'console'
[]
[m11_rate]
type = ParsedPostprocessor
expression = 'fcn*if(p11>0.8,1,if(p11<0.3,0.5,0.2+p11))'
constant_names = 'fcn'
constant_expressions = '8'
pp_names = 'p11'
[]
[m11_expect]
type = ParsedPostprocessor
expression = 'm11_prev-m11_rate*area*dt'
constant_names = 'area dt'
constant_expressions = '0.5 1E-3'
pp_names = 'm11_prev m11_rate'
[]
[]
[BCs]
[flux]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
PT_shift = pt_shift
pt_vals = '0.0 0.5'
multipliers = '0.5 1'
variable = pp
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 8
save_in = flux_out
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1E-3
end_time = 1E-2
nl_rel_tol = 1E-12
nl_abs_tol = 1E-12
[]
[Outputs]
file_base = s04
[console]
type = Console
execute_on = 'nonlinear linear'
[]
[csv]
type = CSV
execute_on = 'timestep_end'
[]
[]
(modules/porous_flow/test/tests/actions/basicthm_thm.i)
# PorousFlowBasicTHM action with coupling_type = ThermoHydroMechanical
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
displacements = 'disp_x disp_y'
biot_coefficient = 1.0
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[temperature]
initial_condition = 293
scaling = 1e-6
[]
[disp_x]
scaling = 1e-6
[]
[disp_y]
scaling = 1e-6
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydroMechanical
gravity = '0 0 0'
fp = simple_fluid
eigenstrain_names = thermal_contribution
use_displaced_mesh = false
add_stress_aux = false
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 313
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[constant_outflow_temperature]
type = DirichletBC
variable = temperature
value = 293
boundary = outflow_area
[]
[top_bottom]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'top bottom'
[]
[right]
type = DirichletBC
variable = disp_x
value = 0
boundary = right
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caprock aquifer'
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5e9
poissons_ratio = 0.0
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = thermal_contribution
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
thermal_expansion_coeff = 0.001
eigenstrain_name = thermal_contribution
stress_free_temperature = 293
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-12
nl_rel_tol = 1E-10
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/actions/basicthm_th.i)
# PorousFlowBasicTHM action with coupling_type = ThermoHydroGenerator
# (no mechanical effects)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[temperature]
initial_condition = 293
scaling = 1e-6
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = simple_fluid
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 313
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[constant_outflow_temperature]
type = DirichletBC
variable = temperature
value = 293
boundary = outflow_area
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
biot_coefficient = 0.8
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caprock aquifer'
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-15
nl_rel_tol = 1e-14
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D_adaptivity.i)
# Using flux-limited TVD advection ala Kuzmin and Turek, mploying PorousFlow Kernels and UserObjects, with superbee flux-limiter
# 1D version with adaptivity
[Mesh]
type = GeneratedMesh
dim = 1
nx = 10
xmin = 0
xmax = 1
[]
[Adaptivity]
initial_steps = 1
initial_marker = tracer_marker
marker = tracer_marker
max_h_level = 1
[Markers]
[tracer_marker]
type = ValueRangeMarker
variable = tracer
lower_bound = 0.02
upper_bound = 0.98
[]
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowFluxLimitedTVDAdvection
variable = tracer
advective_flux_calculator = advective_flux_calculator_0
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowFluxLimitedTVDAdvection
variable = porepressure
advective_flux_calculator = advective_flux_calculator_1
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[advective_flux_calculator_0]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 0
[]
[advective_flux_calculator_1]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 11
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-2
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/examples/restart/gas_injection_new_mesh.i)
# Using the results from the equilibrium run to provide the initial condition for
# porepressure, we now inject a gas phase into the brine-saturated reservoir. In this
# example, the mesh is not identical to the mesh used in gravityeq.i. Rather, it is
# generated so that it is more refined near the injection boundary and at the top of
# the model, as that is where the gas plume will be present.
#
# To use the hydrostatic pressure calculated using the gravity equilibrium run as the initial
# condition for the pressure, a SolutionUserObject is used, along with a SolutionFunction to
# interpolate the pressure from the gravity equilibrium run to the initial condition for liqiud
# porepressure in this example.
#
# Even though the gravity equilibrium is established using a 2D mesh, in this example,
# we use a mesh shifted 0.1 m to the right and rotate it about the Y axis to make a 2D radial
# model.
#
# Methane injection takes place over the surface of the hole created by rotating the mesh,
# and hence the injection area is 2 pi r h. We can calculate this using an AreaPostprocessor,
# and then use this in a ParsedFunction to calculate the injection rate so that 10 kg/s of
# methane is injected.
#
# Note: as this example uses the results from a previous simulation, gravityeq.i MUST be
# run before running this input file.
[Mesh]
type = GeneratedMesh
dim = 2
ny = 25
nx = 50
ymax = 100
xmin = 0.1
xmax = 5000
bias_x = 1.05
bias_y = 0.95
coord_type = RZ
rz_coord_axis = Y
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 -9.81 0'
temperature_unit = Celsius
[]
[Variables]
[pp_liq]
[]
[sat_gas]
initial_condition = 0
[]
[]
[ICs]
[ppliq_ic]
type = FunctionIC
variable = pp_liq
function = ppliq_ic
[]
[]
[AuxVariables]
[temperature]
initial_condition = 50
[]
[xnacl]
initial_condition = 0.1
[]
[brine_density]
family = MONOMIAL
order = CONSTANT
[]
[methane_density]
family = MONOMIAL
order = CONSTANT
[]
[massfrac_ph0_sp0]
initial_condition = 1
[]
[massfrac_ph1_sp0]
initial_condition = 0
[]
[pp_gas]
family = MONOMIAL
order = CONSTANT
[]
[sat_liq]
family = MONOMIAL
order = CONSTANT
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
variable = pp_liq
[]
[flux0]
type = PorousFlowAdvectiveFlux
variable = pp_liq
[]
[mass1]
type = PorousFlowMassTimeDerivative
variable = sat_gas
fluid_component = 1
[]
[flux1]
type = PorousFlowAdvectiveFlux
variable = sat_gas
fluid_component = 1
[]
[]
[AuxKernels]
[brine_density]
type = PorousFlowPropertyAux
property = density
variable = brine_density
execute_on = 'initial timestep_end'
[]
[methane_density]
type = PorousFlowPropertyAux
property = density
variable = methane_density
phase = 1
execute_on = 'initial timestep_end'
[]
[pp_gas]
type = PorousFlowPropertyAux
property = pressure
phase = 1
variable = pp_gas
execute_on = 'initial timestep_end'
[]
[sat_liq]
type = PorousFlowPropertyAux
property = saturation
variable = sat_liq
execute_on = 'initial timestep_end'
[]
[]
[BCs]
[gas_injection]
type = PorousFlowSink
boundary = left
variable = sat_gas
flux_function = injection_rate
fluid_phase = 1
[]
[brine_out]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pp_liq
multipliers = '0 1e9'
pt_vals = '0 1e9'
fluid_phase = 0
flux_function = 1e-6
use_mobility = true
use_relperm = true
mass_fraction_component = 0
[]
[]
[Functions]
[injection_rate]
type = ParsedFunction
symbol_values = injection_area
symbol_names = area
expression = '-1/area'
[]
[ppliq_ic]
type = SolutionFunction
solution = soln
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp_liq sat_gas'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1e-5
m = 0.5
sat_lr = 0.2
pc_max = 1e7
[]
[soln]
type = SolutionUserObject
mesh = gravityeq_out.e
system_variables = porepressure
[]
[]
[FluidProperties]
[brine]
type = BrineFluidProperties
[]
[methane]
type = MethaneFluidProperties
[]
[methane_tab]
type = TabulatedBicubicFluidProperties
fp = methane
save_file = false
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temperature
[]
[ps]
type = PorousFlow2PhasePS
phase0_porepressure = pp_liq
phase1_saturation = sat_gas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
[]
[brine]
type = PorousFlowBrine
compute_enthalpy = false
compute_internal_energy = false
xnacl = xnacl
phase = 0
[]
[methane]
type = PorousFlowSingleComponentFluid
compute_enthalpy = false
compute_internal_energy = false
fp = methane_tab
phase = 1
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 5e-14 0 0 0 1e-13'
[]
[relperm_liq]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.2
sum_s_res = 0.3
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 1
s_res = 0.1
sum_s_res = 0.3
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = ' asm lu NONZERO'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e8
nl_abs_tol = 1e-12
nl_rel_tol = 1e-06
nl_max_its = 20
dtmax = 1e6
[TimeStepper]
type = IterationAdaptiveDT
dt = 1e1
growth_factor = 1.5
[]
[]
[Postprocessors]
[mass_ph0]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'initial timestep_end'
[]
[mass_ph1]
type = PorousFlowFluidMass
fluid_component = 1
execute_on = 'initial timestep_end'
[]
[injection_area]
type = AreaPostprocessor
boundary = left
execute_on = initial
[]
[]
[Outputs]
execute_on = 'initial timestep_end'
exodus = true
perf_graph = true
[]
(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_1D.i)
# Using flux-limited TVD advection ala Kuzmin and Turek, mploying PorousFlow Kernels and UserObjects, with superbee flux-limiter
# 1D version
[Mesh]
type = GeneratedMesh
dim = 1
nx = 10
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowFluxLimitedTVDAdvection
variable = tracer
advective_flux_calculator = advective_flux_calculator_0
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowFluxLimitedTVDAdvection
variable = porepressure
advective_flux_calculator = advective_flux_calculator_1
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[advective_flux_calculator_0]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 0
[]
[advective_flux_calculator_1]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 11
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-2
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/sinks/s09_fully_saturated.i)
# Apply a piecewise-linear sink flux to the right-hand side and watch fluid flow to it
#
# This test has a single phase with two components. The test initialises with
# the porous material fully filled with component=1. The left-hand side is fixed
# at porepressure=1 and mass-fraction of the zeroth component being unity.
# The right-hand side has a very strong piecewise-linear flux that keeps the
# porepressure~0 at that side. Fluid mass is extracted by this flux in proportion
# to the fluid component mass fraction.
#
# Therefore, the zeroth fluid component will flow from left to right (down the
# pressure gradient).
#
# The important DE is
# porosity * dc/dt = (perm / visc) * grad(P) * grad(c)
# which is true for c = mass-fraction, and very large bulk modulus of the fluid.
# For grad(P) constant in time and space (as in this example) this is just the
# advection equation for c, with velocity = perm / visc / porosity. The parameters
# are chosen to velocity = 1 m/s.
# In the numerical world, and especially with full upwinding, the advection equation
# suffers from diffusion. In this example, the diffusion is obvious when plotting
# the mass-fraction along the line, but the average velocity of the front is still
# correct at 1 m/s.
# This test uses the FullySaturated version of the flow Kernel. This does not
# suffer from as much numerical diffusion as the standard PorousFlow Kernel since
# it does not employ any upwinding.
[Mesh]
type = GeneratedMesh
dim = 1
nx = 100
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp frac'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[Variables]
[pp]
[]
[frac]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
function = 1-x
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = frac
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pp
[]
[flux0]
type = PorousFlowFullySaturatedDarcyFlow
fluid_component = 0
gravity = '0 0 0'
variable = frac
[]
[flux1]
type = PorousFlowFullySaturatedDarcyFlow
fluid_component = 1
gravity = '0 0 0'
variable = pp
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1e10 # need large in order for constant-velocity advection
density0 = 1 # almost irrelevant, except that the ability of the right BC to keep P fixed at zero is related to density_P0
thermal_expansion = 0
viscosity = 11
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = frac
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.1 0 0 0 1.1 0 0 0 1.1'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 2 # irrelevant in this fully-saturated situation
phase = 0
[]
[]
[BCs]
[lhs_fixed_a]
type = DirichletBC
boundary = 'left'
variable = frac
value = 1
[]
[lhs_fixed_b]
type = DirichletBC
boundary = 'left'
variable = pp
value = 1
[]
[flux0]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '-100 100'
multipliers = '-1 1'
variable = frac # the zeroth comonent
mass_fraction_component = 0
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 1E4
[]
[flux1]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
pt_vals = '-100 100'
multipliers = '-1 1'
variable = pp # comonent 1
mass_fraction_component = 1
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 1E4
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu 10000 NONZERO 2'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1E-2
end_time = 1
nl_rel_tol = 1E-11
nl_abs_tol = 1E-11
[]
[VectorPostprocessors]
[mf]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 100
sort_by = x
variable = frac
[]
[]
[Outputs]
file_base = s09_fully_saturated
[console]
type = Console
execute_on = 'nonlinear linear'
[]
[csv]
type = CSV
sync_times = '0.1 0.5 1'
sync_only = true
[]
time_step_interval = 10
[]
(modules/porous_flow/examples/restart/gas_injection.i)
# Using the results from the equilibrium run to provide the initial condition for
# porepressure, we now inject a gas phase into the brine-saturated reservoir. In this
# example, where the mesh used is identical to the mesh used in gravityeq.i, we can use
# the basic restart capability by simply setting the initial condition for porepressure
# using the results from gravityeq.i.
#
# Even though the gravity equilibrium is established using a 2D mesh, in this example,
# we shift the mesh 0.1 m to the right and rotate it about the Y axis to make a 2D radial
# model.
#
# Methane injection takes place over the surface of the hole created by rotating the mesh,
# and hence the injection area is 2 pi r h. We can calculate this using an AreaPostprocessor,
# and then use this in a ParsedFunction to calculate the injection rate so that 10 kg/s of
# methane is injected.
#
# Results can be improved by uniformly refining the initial mesh.
#
# Note: as this example uses the results from a previous simulation, gravityeq.i MUST be
# run before running this input file.
[Mesh]
uniform_refine = 1
[file]
type = FileMeshGenerator
file = gravityeq_out.e
[]
[translate]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0.1 0 0'
input = file
[]
coord_type = RZ
rz_coord_axis = Y
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 -9.81 0'
temperature_unit = Celsius
[]
[Variables]
[pp_liq]
initial_from_file_var = porepressure
[]
[sat_gas]
initial_condition = 0
[]
[]
[AuxVariables]
[temperature]
initial_condition = 50
[]
[xnacl]
initial_condition = 0.1
[]
[brine_density]
family = MONOMIAL
order = CONSTANT
[]
[methane_density]
family = MONOMIAL
order = CONSTANT
[]
[massfrac_ph0_sp0]
initial_condition = 1
[]
[massfrac_ph1_sp0]
initial_condition = 0
[]
[pp_gas]
family = MONOMIAL
order = CONSTANT
[]
[sat_liq]
family = MONOMIAL
order = CONSTANT
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
variable = pp_liq
[]
[flux0]
type = PorousFlowAdvectiveFlux
variable = pp_liq
[]
[mass1]
type = PorousFlowMassTimeDerivative
variable = sat_gas
fluid_component = 1
[]
[flux1]
type = PorousFlowAdvectiveFlux
variable = sat_gas
fluid_component = 1
[]
[]
[AuxKernels]
[brine_density]
type = PorousFlowPropertyAux
property = density
variable = brine_density
execute_on = 'initial timestep_end'
[]
[methane_density]
type = PorousFlowPropertyAux
property = density
variable = methane_density
phase = 1
execute_on = 'initial timestep_end'
[]
[pp_gas]
type = PorousFlowPropertyAux
property = pressure
phase = 1
variable = pp_gas
execute_on = 'initial timestep_end'
[]
[sat_liq]
type = PorousFlowPropertyAux
property = saturation
variable = sat_liq
execute_on = 'initial timestep_end'
[]
[]
[BCs]
[gas_injection]
type = PorousFlowSink
boundary = left
variable = sat_gas
flux_function = injection_rate
fluid_phase = 1
[]
[brine_out]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pp_liq
multipliers = '0 1e9'
pt_vals = '0 1e9'
fluid_phase = 0
flux_function = 1e-6
use_mobility = true
[]
[]
[Functions]
[injection_rate]
type = ParsedFunction
symbol_values = injection_area
symbol_names = area
expression = '-10/area'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp_liq sat_gas'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1e-5
m = 0.5
sat_lr = 0.2
[]
[]
[FluidProperties]
[brine]
type = BrineFluidProperties
[]
[methane]
type = MethaneFluidProperties
[]
[methane_tab]
type = TabulatedBicubicFluidProperties
fp = methane
save_file = false
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temperature
[]
[ps]
type = PorousFlow2PhasePS
phase0_porepressure = pp_liq
phase1_saturation = sat_gas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'massfrac_ph0_sp0 massfrac_ph1_sp0'
[]
[brine]
type = PorousFlowBrine
compute_enthalpy = false
compute_internal_energy = false
xnacl = xnacl
phase = 0
[]
[methane]
type = PorousFlowSingleComponentFluid
compute_enthalpy = false
compute_internal_energy = false
fp = methane_tab
phase = 1
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[relperm_liq]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.2
sum_s_res = 0.3
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 1
s_res = 0.1
sum_s_res = 0.3
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = ' asm lu NONZERO'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e8
nl_abs_tol = 1e-12
nl_rel_tol = 1e-06
nl_max_its = 20
dtmax = 1e6
[TimeStepper]
type = IterationAdaptiveDT
dt = 1e1
[]
[]
[Postprocessors]
[mass_ph0]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'initial timestep_end'
[]
[mass_ph1]
type = PorousFlowFluidMass
fluid_component = 1
execute_on = 'initial timestep_end'
[]
[injection_area]
type = AreaPostprocessor
boundary = left
execute_on = initial
[]
[]
[Outputs]
execute_on = 'initial timestep_end'
exodus = true
perf_graph = true
checkpoint = true
[]
(modules/porous_flow/test/tests/sinks/PorousFlowPiecewiseLinearSink_BC_eg1.i)
## This is an example input file showing how to set a Type I (Dirichlet) BC with PorousFlowPiecewiseLinearSink
##
## Problem setup:
## - The boundaries are set to P(x = 0) = 2e6 Pa, P(x = 1) = 1e6 and run to steady state.
## - The 2d domain is 1 m x 1 m
## - The permeability is set to 1E-15 m2, fluid viscosity = 1E-3 Pa-s
## - The steady state flux is calculated q = -k/mu*grad(P) = 1e-6 m/s
##
## Problem verification (in csv output):
## - The flux in and out of the domain are 1e-6 m/s (matching steady state solution)
## - The pressure at the left and right boundaries are set to 2e6 and 1e6 Pa, respectively
[Mesh]
type = GeneratedMesh
dim = 2
nx = 5
xmin = 0
xmax = 1
ny = 2
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1.5e6 # initial pressure in domain
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
coupling_type = Hydro
gravity = '0 0 0'
fp = the_simple_fluid
[]
[AuxVariables]
[fluxes_out]
[]
[fluxes_in]
[]
[]
[BCs]
[in_left]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = 'left'
pt_vals = '-1e9 1e9' # x coordinates defining g
multipliers = '-1e9 1e9' # y coordinates defining g
PT_shift = 2.E6 # BC pressure
flux_function = 1E-5 # Variable C
fluid_phase = 0
save_in = fluxes_out
[]
[out_right]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = 'right'
pt_vals = '-1e9 1e9' # x coordinates defining g
multipliers = '-1e9 1e9' # y coordinates defining g
PT_shift = 1.E6 # BC pressure
flux_function = 1E-6 # Variable C
fluid_phase = 0
save_in = fluxes_in
[]
[]
[Postprocessors]
[left_flux]
type = NodalSum
boundary = 'left'
variable = fluxes_out
execute_on = 'timestep_end'
[]
[right_flux]
type = NodalSum
boundary = 'right'
variable = fluxes_in
execute_on = 'timestep_end'
[]
[left_pressure]
type = SideAverageValue
boundary = 'left'
variable = porepressure
execute_on = 'timestep_end'
[]
[right_pressure]
type = SideAverageValue
boundary = 'right'
variable = porepressure
execute_on = 'timestep_end'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0
viscosity = 1.0E-3
density0 = 1000.0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-10
[]
[Outputs]
csv = true
[]
(modules/porous_flow/examples/co2_intercomparison/1Dradial/1Dradial.i)
# Intercomparison problem 3: Radial flow from an injection well
#
# From Pruess et al, Code intercomparison builds confidence in
# numerical simulation models for geologic disposal of CO2, Energy 29 (2004)
#
# A variation with zero salinity can be run by changing the initial condition
# of the AuxVariable xnacl
[Mesh]
type = GeneratedMesh
dim = 1
nx = 500
xmax = 10000
bias_x = 1.01
coord_type = 'RZ'
rz_coord_axis = Y
[]
[GlobalParams]
PorousFlowDictator = 'dictator'
gravity = '0 0 0'
[]
[AuxVariables]
[pressure_liquid]
order = CONSTANT
family = MONOMIAL
[]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[x1]
order = CONSTANT
family = MONOMIAL
[]
[y0]
order = CONSTANT
family = MONOMIAL
[]
[xnacl]
initial_condition = 0.15
[]
[]
[AuxKernels]
[pressure_liquid]
type = PorousFlowPropertyAux
variable = pressure_liquid
property = pressure
phase = 0
execute_on = 'timestep_end'
[]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = 'timestep_end'
[]
[x1]
type = PorousFlowPropertyAux
variable = x1
property = mass_fraction
phase = 0
fluid_component = 1
execute_on = 'timestep_end'
[]
[y0]
type = PorousFlowPropertyAux
variable = y0
property = mass_fraction
phase = 1
fluid_component = 0
execute_on = 'timestep_end'
[]
[]
[Variables]
[pgas]
initial_condition = 12e6
[]
[zi]
initial_condition = 0
scaling = 1e4
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pgas
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pgas
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = zi
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = zi
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi'
number_fluid_phases = 2
number_fluid_components = 3
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 5.099e-5
m = 0.457
sat_lr = 0.0
pc_max = 1e7
[]
[fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2
capillary_pressure = pc
[]
[]
[FluidProperties]
[co2sw]
type = CO2FluidProperties
[]
[co2]
type = TabulatedBicubicFluidProperties
fp = co2sw
[]
[water]
type = Water97FluidProperties
[]
[watertab]
type = TabulatedBicubicFluidProperties
fp = water
temperature_min = 273.15
temperature_max = 573.15
fluid_property_output_file = water_fluid_properties.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water_fluid_properties.csv
[]
[brine]
type = BrineFluidProperties
water_fp = watertab
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = '45'
[]
[brineco2]
type = PorousFlowFluidState
gas_porepressure = 'pgas'
z = 'zi'
temperature_unit = Celsius
xnacl = 'xnacl'
capillary_pressure = pc
fluid_state = fs
[]
[porosity]
type = PorousFlowPorosityConst
porosity = '0.12'
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityVG
m = 0.457
phase = 0
s_res = 0.3
sum_s_res = 0.35
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 1
s_res = 0.05
sum_s_res = 0.35
[]
[]
[BCs]
[rightwater]
type = PorousFlowPiecewiseLinearSink
boundary = 'right'
variable = pgas
use_mobility = true
PorousFlowDictator = dictator
fluid_phase = 0
multipliers = '0 1e9'
PT_shift = '12e6'
pt_vals = '0 1e9'
mass_fraction_component = 0
use_relperm = true
[]
[rightco2]
type = PorousFlowPiecewiseLinearSink
variable = zi
boundary = 'right'
use_mobility = true
PorousFlowDictator = dictator
fluid_phase = 1
multipliers = '0 1e9'
PT_shift = '12e6'
pt_vals = '0 1e9'
mass_fraction_component = 1
use_relperm = true
[]
[]
[DiracKernels]
[source]
type = PorousFlowSquarePulsePointSource
point = '0 0 0'
mass_flux = 1
variable = zi
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres bjacobi lu NONZERO'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 8.64e8
nl_max_its = 25
l_max_its = 100
dtmax = 5e6
[TimeStepper]
type = IterationAdaptiveDT
dt = 100
[]
[]
[VectorPostprocessors]
[vars]
type = NodalValueSampler
sort_by = x
variable = 'pgas zi xnacl'
execute_on = 'timestep_end'
outputs = spatial
[]
[auxvars]
type = ElementValueSampler
sort_by = x
variable = 'saturation_gas x1 y0'
execute_on = 'timestep_end'
outputs = spatial
[]
[]
[Postprocessors]
[pgas]
type = PointValue
point = '25.25 0 0'
variable = pgas
outputs = time
[]
[sgas]
type = PointValue
point = '25.25 0 0'
variable = saturation_gas
outputs = time
[]
[zi]
type = PointValue
point = '25.25 0 0'
variable = zi
outputs = time
[]
[massgas]
type = PorousFlowFluidMass
fluid_component = 1
outputs = time
[]
[x1]
type = PointValue
point = '25.25 0 0'
variable = x1
outputs = time
[]
[y0]
type = PointValue
point = '25.25 0 0'
variable = y0
outputs = time
[]
[xnacl]
type = PointValue
point = '25.25 0 0'
variable = xnacl
outputs = time
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
sync_times = '2.592e6 8.64e6 8.64e7 8.64e8'
[time]
type = CSV
[]
[spatial]
type = CSV
sync_only = true
[]
[]
(modules/porous_flow/examples/tutorial/11.i)
# Two-phase borehole injection problem
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
dmin = 0
dmax = 90
[]
[make3D]
input = annular
type = MeshExtruderGenerator
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
[]
[shift_down]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0 0 -6'
input = make3D
[]
[aquifer]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
input = shift_down
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pwater pgas T disp_x disp_y'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
gravity = '0 0 0'
biot_coefficient = 1.0
PorousFlowDictator = dictator
[]
[Variables]
[pwater]
initial_condition = 20E6
[]
[pgas]
initial_condition = 20.1E6
[]
[T]
initial_condition = 330
scaling = 1E-5
[]
[disp_x]
scaling = 1E-5
[]
[disp_y]
scaling = 1E-5
[]
[]
[Kernels]
[mass_water_dot]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pwater
[]
[flux_water]
type = PorousFlowAdvectiveFlux
fluid_component = 0
use_displaced_mesh = false
variable = pwater
[]
[vol_strain_rate_water]
type = PorousFlowMassVolumetricExpansion
fluid_component = 0
variable = pwater
[]
[mass_co2_dot]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pgas
[]
[flux_co2]
type = PorousFlowAdvectiveFlux
fluid_component = 1
use_displaced_mesh = false
variable = pgas
[]
[vol_strain_rate_co2]
type = PorousFlowMassVolumetricExpansion
fluid_component = 1
variable = pgas
[]
[energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = T
[]
[advection]
type = PorousFlowHeatAdvection
use_displaced_mesh = false
variable = T
[]
[conduction]
type = PorousFlowHeatConduction
use_displaced_mesh = false
variable = T
[]
[vol_strain_rate_heat]
type = PorousFlowHeatVolumetricExpansion
variable = T
[]
[grad_stress_x]
type = StressDivergenceTensors
temperature = T
variable = disp_x
eigenstrain_names = thermal_contribution
use_displaced_mesh = false
component = 0
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
variable = disp_x
use_displaced_mesh = false
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
temperature = T
variable = disp_y
eigenstrain_names = thermal_contribution
use_displaced_mesh = false
component = 1
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
variable = disp_y
use_displaced_mesh = false
component = 1
[]
[]
[AuxVariables]
[disp_z]
[]
[effective_fluid_pressure]
family = MONOMIAL
order = CONSTANT
[]
[mass_frac_phase0_species0]
initial_condition = 1 # all water in phase=0
[]
[mass_frac_phase1_species0]
initial_condition = 0 # no water in phase=1
[]
[sgas]
family = MONOMIAL
order = CONSTANT
[]
[swater]
family = MONOMIAL
order = CONSTANT
[]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_tt]
family = MONOMIAL
order = CONSTANT
[]
[stress_zz]
family = MONOMIAL
order = CONSTANT
[]
[porosity]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[effective_fluid_pressure]
type = ParsedAux
coupled_variables = 'pwater pgas swater sgas'
expression = 'pwater * swater + pgas * sgas'
variable = effective_fluid_pressure
[]
[swater]
type = PorousFlowPropertyAux
variable = swater
property = saturation
phase = 0
execute_on = timestep_end
[]
[sgas]
type = PorousFlowPropertyAux
variable = sgas
property = saturation
phase = 1
execute_on = timestep_end
[]
[stress_rr]
type = RankTwoScalarAux
variable = stress_rr
rank_two_tensor = stress
scalar_type = RadialStress
point1 = '0 0 0'
point2 = '0 0 1'
execute_on = timestep_end
[]
[stress_tt]
type = RankTwoScalarAux
variable = stress_tt
rank_two_tensor = stress
scalar_type = HoopStress
point1 = '0 0 0'
point2 = '0 0 1'
execute_on = timestep_end
[]
[stress_zz]
type = RankTwoAux
variable = stress_zz
rank_two_tensor = stress
index_i = 2
index_j = 2
execute_on = timestep_end
[]
[porosity]
type = PorousFlowPropertyAux
variable = porosity
property = porosity
execute_on = timestep_end
[]
[]
[BCs]
[roller_tmax]
type = DirichletBC
variable = disp_x
value = 0
boundary = dmax
[]
[roller_tmin]
type = DirichletBC
variable = disp_y
value = 0
boundary = dmin
[]
[pinned_top_bottom_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'top bottom'
[]
[pinned_top_bottom_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'top bottom'
[]
[cavity_pressure_x]
type = Pressure
boundary = injection_area
variable = disp_x
component = 0
postprocessor = constrained_effective_fluid_pressure_at_wellbore
use_displaced_mesh = false
[]
[cavity_pressure_y]
type = Pressure
boundary = injection_area
variable = disp_y
component = 1
postprocessor = constrained_effective_fluid_pressure_at_wellbore
use_displaced_mesh = false
[]
[cold_co2]
type = DirichletBC
boundary = injection_area
variable = T
value = 290 # injection temperature
use_displaced_mesh = false
[]
[constant_co2_injection]
type = PorousFlowSink
boundary = injection_area
variable = pgas
fluid_phase = 1
flux_function = -1E-4
use_displaced_mesh = false
[]
[outer_water_removal]
type = PorousFlowPiecewiseLinearSink
boundary = rmax
variable = pwater
fluid_phase = 0
pt_vals = '0 1E9'
multipliers = '0 1E8'
PT_shift = 20E6
use_mobility = true
use_relperm = true
use_displaced_mesh = false
[]
[outer_co2_removal]
type = PorousFlowPiecewiseLinearSink
boundary = rmax
variable = pgas
fluid_phase = 1
pt_vals = '0 1E9'
multipliers = '0 1E8'
PT_shift = 20.1E6
use_mobility = true
use_relperm = true
use_displaced_mesh = false
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[tabulated_water]
type = TabulatedFluidProperties
fp = true_water
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = water97_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated_11.csv
[]
[true_co2]
type = CO2FluidProperties
[]
[tabulated_co2]
type = TabulatedFluidProperties
fp = true_co2
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = co2_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = co2_tabulated_11.csv
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = T
[]
[saturation_calculator]
type = PorousFlow2PhasePP
phase0_porepressure = pwater
phase1_porepressure = pgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
[]
[water]
type = PorousFlowSingleComponentFluid
fp = tabulated_water
phase = 0
[]
[co2]
type = PorousFlowSingleComponentFluid
fp = tabulated_co2
phase = 1
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 4
s_res = 0.1
sum_s_res = 0.2
phase = 0
[]
[relperm_co2]
type = PorousFlowRelativePermeabilityBC
nw_phase = true
lambda = 2
s_res = 0.1
sum_s_res = 0.2
phase = 1
[]
[porosity_mat]
type = PorousFlowPorosity
fluid = true
mechanical = true
thermal = true
porosity_zero = 0.1
reference_temperature = 330
reference_porepressure = 20E6
thermal_expansion_coeff = 15E-6 # volumetric
solid_bulk = 8E9 # unimportant since biot = 1
[]
[permeability_aquifer]
type = PorousFlowPermeabilityKozenyCarman
block = aquifer
poroperm_function = kozeny_carman_phi0
phi0 = 0.1
n = 2
m = 2
k0 = 1E-12
[]
[permeability_caps]
type = PorousFlowPermeabilityKozenyCarman
block = caps
poroperm_function = kozeny_carman_phi0
phi0 = 0.1
n = 2
m = 2
k0 = 1E-15
k_anisotropy = '1 0 0 0 1 0 0 0 0.1'
[]
[rock_thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2 0 0 0 2 0 0 0 2'
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 1100
density = 2300
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5E9
poissons_ratio = 0.0
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = 'thermal_contribution initial_stress'
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = T
thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
eigenstrain_name = thermal_contribution
stress_free_temperature = 330
[]
[initial_strain]
type = ComputeEigenstrainFromInitialStress
initial_stress = '20E6 0 0 0 20E6 0 0 0 20E6'
eigenstrain_name = initial_stress
[]
[stress]
type = ComputeLinearElasticStress
[]
[effective_fluid_pressure_mat]
type = PorousFlowEffectiveFluidPressure
[]
[volumetric_strain]
type = PorousFlowVolumetricStrain
[]
[]
[Postprocessors]
[effective_fluid_pressure_at_wellbore]
type = PointValue
variable = effective_fluid_pressure
point = '1 0 0'
execute_on = timestep_begin
use_displaced_mesh = false
[]
[constrained_effective_fluid_pressure_at_wellbore]
type = FunctionValuePostprocessor
function = constrain_effective_fluid_pressure
execute_on = timestep_begin
[]
[]
[Functions]
[constrain_effective_fluid_pressure]
type = ParsedFunction
symbol_names = effective_fluid_pressure_at_wellbore
symbol_values = effective_fluid_pressure_at_wellbore
expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E3
[TimeStepper]
type = IterationAdaptiveDT
dt = 1E3
growth_factor = 1.2
optimal_iterations = 10
[]
nl_abs_tol = 1E-7
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/newton_cooling/nc02.i)
# Newton cooling from a bar. 1-phase steady
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1000
ny = 1
xmin = 0
xmax = 100
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pressure'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Variables]
[pressure]
[]
[]
[ICs]
[pressure]
type = FunctionIC
variable = pressure
function = '(2-x/100)*1E6'
[]
[]
[Kernels]
[flux]
type = PorousFlowAdvectiveFlux
fluid_component = 0
gravity = '0 0 0'
variable = pressure
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1e6
density0 = 1000
thermal_expansion = 0
viscosity = 1e-3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
n = 2
phase = 0
[]
[]
[BCs]
[left]
type = DirichletBC
variable = pressure
boundary = left
value = 2E6
[]
[newton]
type = PorousFlowPiecewiseLinearSink
variable = pressure
boundary = right
pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 1
[]
[]
[VectorPostprocessors]
[porepressure]
type = LineValueSampler
variable = pressure
start_point = '0 0.5 0'
end_point = '100 0.5 0'
sort_by = x
num_points = 20
execute_on = timestep_end
[]
[]
[Preconditioning]
active = 'andy'
[andy]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
petsc_options_value = 'gmres asm lu 100 NONZERO 2 1E-12 1E-15'
[]
[]
[Executioner]
type = Steady
[]
[Outputs]
file_base = nc02
execute_on = timestep_end
exodus = false
[along_line]
type = CSV
execute_vector_postprocessors_on = timestep_end
[]
[]
(modules/porous_flow/test/tests/recover/pffltvd.i)
# Tests that PorousFlow can successfully recover using a checkpoint file.
# This test contains stateful material properties, adaptivity, integrated
# boundary conditions with nodal-sized materials, and TVD flux limiting.
#
# This test file is run three times:
# 1) The full input file is run to completion
# 2) The input file is run for half the time and checkpointing is included
# 3) The input file is run in recovery using the checkpoint data
#
# The final output of test 3 is compared to the final output of test 1 to verify
# that recovery was successful.
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
nx = 10
xmin = 0
xmax = 1
[]
# To get consistent ordering of results with distributed meshes
allow_renumbering = false
[]
[Adaptivity]
initial_steps = 1
initial_marker = tracer_marker
marker = tracer_marker
max_h_level = 1
[Markers]
[tracer_marker]
type = ValueRangeMarker
variable = tracer
lower_bound = 0.02
upper_bound = 0.98
[]
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '2 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowFluxLimitedTVDAdvection
variable = tracer
advective_flux_calculator = advective_flux_calculator_0
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowFluxLimitedTVDAdvection
variable = porepressure
advective_flux_calculator = advective_flux_calculator_1
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 2
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[advective_flux_calculator_0]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 0
[]
[advective_flux_calculator_1]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
[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'
[]
[]
[VectorPostprocessors]
[tracer]
type = NodalValueSampler
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 0.2
dt = 0.05
[]
[Outputs]
csv = true
[]
(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_3D.i)
# Using flux-limited TVD advection ala Kuzmin and Turek, employing PorousFlow Kernels and UserObjects, with superbee flux-limiter
# 3D version
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
xmin = 0
xmax = 1
ny = 4
ymin = 0
ymax = 0.5
nz = 3
zmin = 0
zmax = 2
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowFluxLimitedTVDAdvection
variable = tracer
advective_flux_calculator = advective_flux_calculator_0
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowFluxLimitedTVDAdvection
variable = porepressure
advective_flux_calculator = advective_flux_calculator_1
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[advective_flux_calculator_0]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 0
[]
[advective_flux_calculator_1]
type = PorousFlowAdvectiveFluxCalculatorSaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0.5 2'
num_points = 11
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 0.3
dt = 6E-2
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
print_linear_residuals = false
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/numerical_diffusion/no_action.i)
# Using upwinded and mass-lumped PorousFlow Kernels: this is equivalent of fully_saturated_action.i with stabilization = Full
[Mesh]
type = GeneratedMesh
dim = 1
nx = 100
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = tracer
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = porepressure
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 101
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-1
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/newton_cooling/nc04.i)
# Newton cooling from a bar. Heat conduction
[Mesh]
type = GeneratedMesh
dim = 2
nx = 100
ny = 1
xmin = 0
xmax = 100
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp'
number_fluid_phases = 0
number_fluid_components = 0
[]
[]
[Variables]
[temp]
[]
[]
[ICs]
[temp]
type = FunctionIC
variable = temp
function = '2-x/100'
[]
[]
[Kernels]
[conduction]
type = PorousFlowHeatConduction
variable = temp
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[thermal_conductivity_irrelevant]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1E2 0 0 0 1E2 0 0 0 1E2'
[]
[]
[BCs]
[left]
type = DirichletBC
variable = temp
boundary = left
value = 2
[]
[newton]
type = PorousFlowPiecewiseLinearSink
variable = temp
boundary = right
pt_vals = '0 1 2'
multipliers = '-1 0 1'
flux_function = 1
[]
[]
[VectorPostprocessors]
[temp]
type = LineValueSampler
variable = temp
start_point = '0 0.5 0'
end_point = '100 0.5 0'
sort_by = x
num_points = 11
execute_on = timestep_end
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
petsc_options_value = 'gmres asm lu 100 NONZERO 2 1E-14 1E-12'
[]
[]
[Executioner]
type = Steady
[]
[Outputs]
file_base = nc04
execute_on = timestep_end
exodus = false
[along_line]
type = CSV
execute_vector_postprocessors_on = timestep_end
[]
[]
(modules/porous_flow/test/tests/numerical_diffusion/pffltvd_action.i)
# Using flux-limited TVD advection ala Kuzmin and Turek, employing PorousFlow Kernels and UserObjects, with superbee flux-limiter
# Using the PorousFlowFullySaturated Action
[Mesh]
type = GeneratedMesh
dim = 1
nx = 100
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[PorousFlowFullySaturated]
porepressure = porepressure
coupling_type = Hydro
fp = the_simple_fluid
mass_fraction_vars = tracer
stabilization = KT
flux_limiter_type = superbee
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 101
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-2
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
file_base = pffltvd_out
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/flux_limited_TVD_pflow/pffltvd_2D_trimesh.i)
# Using flux-limited TVD advection ala Kuzmin and Turek, mploying PorousFlow Kernels and UserObjects, with superbee flux-limiter
# 2D version
[Mesh]
type = FileMesh
file = trimesh.msh
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
block = '50'
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.305,0,1))'
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = tracer
[]
[flux0]
type = PorousFlowFluxLimitedTVDAdvection
variable = tracer
advective_flux_calculator = advective_flux_calculator_0
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = porepressure
[]
[flux1]
type = PorousFlowFluxLimitedTVDAdvection
variable = porepressure
advective_flux_calculator = advective_flux_calculator_1
[]
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure tracer'
number_fluid_phases = 1
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[advective_flux_calculator_0]
type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 0
[]
[advective_flux_calculator_1]
type = PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent
flux_limiter_type = superbee
fluid_component = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = tracer
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[relperm]
type = PorousFlowRelativePermeabilityConst
phase = 0
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0.04 0'
num_points = 101
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-2
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
print_linear_residuals = false
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/examples/coal_mining/fine_with_fluid.i)
#################################################################
#
# NOTE:
# The mesh for this model is too large for the MOOSE repository
# so is kept in the the large_media submodule
#
#################################################################
#
# Strata deformation and fluid flow aaround a coal mine - 3D model
#
# A "half model" is used. The mine is 400m deep and
# just the roof is studied (-400<=z<=0). The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long. The outer boundaries
# are 1km from the excavation boundaries.
#
# The excavation takes 0.5 years.
#
# The boundary conditions for this simulation are:
# - disp_x = 0 at x=0 and x=1150
# - disp_y = 0 at y=-1000 and y=1000
# - disp_z = 0 at z=-400, but there is a time-dependent
# Young modulus that simulates excavation
# - wc_x = 0 at y=-1000 and y=1000
# - wc_y = 0 at x=0 and x=1150
# - no flow at x=0, z=-400 and z=0
# - fixed porepressure at y=-1000, y=1000 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# A single-phase unsaturated fluid is used.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa, and time units are measured in years.
#
# The initial porepressure is hydrostatic with P=0 at z=0, so
# Porepressure ~ - 0.01*z MPa, where the fluid has density 1E3 kg/m^3 and
# gravity = = 10 m.s^-2 = 1E-5 MPa m^2/kg.
# To be more accurate, i use
# Porepressure = -bulk * log(1 + g*rho0*z/bulk)
# where bulk=2E3 MPa and rho0=1Ee kg/m^3.
# The initial stress is consistent with the weight force from undrained
# density 2500 kg/m^3, and fluid porepressure, and a Biot coefficient of 0.7, ie,
# stress_zz^effective = 0.025*z + 0.7 * initial_porepressure
# The maximum and minimum principal horizontal effective stresses are
# assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 2 MPa
# MC friction angle = 35 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
# Fluid density at zero porepressure = 1E3 kg/m^3
# Fluid bulk modulus = 2E3 MPa
# Fluid viscosity = 1.1E-3 Pa.s = 1.1E-9 MPa.s = 3.5E-17 MPa.year
#
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
PorousFlowDictator = dictator
biot_coefficient = 0.7
[]
[Mesh]
[file]
type = FileMeshGenerator
file = fine.e
[]
[xmin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = xmin
normal = '-1 0 0'
input = file
[]
[xmax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = xmax
normal = '1 0 0'
input = xmin
[]
[ymin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = ymin
normal = '0 -1 0'
input = xmax
[]
[ymax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = ymax
normal = '0 1 0'
input = ymin
[]
[zmax]
type = SideSetsAroundSubdomainGenerator
block = 30
new_boundary = zmax
normal = '0 0 1'
input = ymax
[]
[zmin]
type = SideSetsAroundSubdomainGenerator
block = 2
new_boundary = zmin
normal = '0 0 -1'
input = zmax
[]
[excav]
type = SubdomainBoundingBoxGenerator
input = zmin
block_id = 1
bottom_left = '0 0 -400'
top_right = '150 1000 -397'
[]
[roof]
type = SideSetsBetweenSubdomainsGenerator
primary_block = 3
paired_block = 1
input = excav
new_boundary = roof
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[wc_x]
[]
[wc_y]
[]
[porepressure]
scaling = 1E-5
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = ini_pp
[]
[]
[Kernels]
[cx_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_x
component = 0
[]
[cy_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_y
component = 1
[]
[cz_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_z
component = 2
[]
[x_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_x
displacements = 'wc_x wc_y wc_z'
component = 0
base_name = couple
[]
[y_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_y
displacements = 'wc_x wc_y wc_z'
component = 1
base_name = couple
[]
[x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[]
[y_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_y
component = 1
[]
[gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
component = 2
variable = disp_z
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
use_displaced_mesh = false
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
variable = porepressure
fluid_component = 0
[]
[mass0]
type = PorousFlowMassTimeDerivative
use_displaced_mesh = false
fluid_component = 0
variable = porepressure
[]
[flux]
type = PorousFlowAdvectiveFlux
use_displaced_mesh = false
variable = porepressure
gravity = '0 0 -10E-6'
fluid_component = 0
[]
[]
[AuxVariables]
[saturation]
order = CONSTANT
family = MONOMIAL
[]
[darcy_x]
order = CONSTANT
family = MONOMIAL
[]
[darcy_y]
order = CONSTANT
family = MONOMIAL
[]
[darcy_z]
order = CONSTANT
family = MONOMIAL
[]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[wc_z]
[]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[stress_yx]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_zx]
order = CONSTANT
family = MONOMIAL
[]
[stress_zy]
order = CONSTANT
family = MONOMIAL
[]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_xx]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_xy]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_xz]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_yx]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_yz]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_zx]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_zy]
order = CONSTANT
family = MONOMIAL
[]
[total_strain_zz]
order = CONSTANT
family = MONOMIAL
[]
[perm_xx]
order = CONSTANT
family = MONOMIAL
[]
[perm_yy]
order = CONSTANT
family = MONOMIAL
[]
[perm_zz]
order = CONSTANT
family = MONOMIAL
[]
[mc_shear]
order = CONSTANT
family = MONOMIAL
[]
[mc_tensile]
order = CONSTANT
family = MONOMIAL
[]
[wp_shear]
order = CONSTANT
family = MONOMIAL
[]
[wp_tensile]
order = CONSTANT
family = MONOMIAL
[]
[wp_shear_f]
order = CONSTANT
family = MONOMIAL
[]
[wp_tensile_f]
order = CONSTANT
family = MONOMIAL
[]
[mc_shear_f]
order = CONSTANT
family = MONOMIAL
[]
[mc_tensile_f]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[saturation_water]
type = PorousFlowPropertyAux
variable = saturation
property = saturation
phase = 0
execute_on = timestep_end
[]
[darcy_x]
type = PorousFlowDarcyVelocityComponent
variable = darcy_x
gravity = '0 0 -10E-6'
component = x
[]
[darcy_y]
type = PorousFlowDarcyVelocityComponent
variable = darcy_y
gravity = '0 0 -10E-6'
component = y
[]
[darcy_z]
type = PorousFlowDarcyVelocityComponent
variable = darcy_z
gravity = '0 0 -10E-6'
component = z
[]
[porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
execute_on = timestep_end
[]
[stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[]
[stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[]
[stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[]
[stress_yx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[]
[stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[]
[stress_zx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[]
[stress_zy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[]
[total_strain_xx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[]
[total_strain_xy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[]
[total_strain_xz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[]
[total_strain_yx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[]
[total_strain_yy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[]
[total_strain_yz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[]
[total_strain_zx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[]
[total_strain_zy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[]
[total_strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[]
[perm_xx]
type = PorousFlowPropertyAux
property = permeability
variable = perm_xx
row = 0
column = 0
execute_on = timestep_end
[]
[perm_yy]
type = PorousFlowPropertyAux
property = permeability
variable = perm_yy
row = 1
column = 1
execute_on = timestep_end
[]
[perm_zz]
type = PorousFlowPropertyAux
property = permeability
variable = perm_zz
row = 2
column = 2
execute_on = timestep_end
[]
[mc_shear]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_internal_parameter
variable = mc_shear
execute_on = timestep_end
[]
[mc_tensile]
type = MaterialStdVectorAux
index = 1
property = mc_plastic_internal_parameter
variable = mc_tensile
execute_on = timestep_end
[]
[wp_shear]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_internal_parameter
variable = wp_shear
execute_on = timestep_end
[]
[wp_tensile]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_internal_parameter
variable = wp_tensile
execute_on = timestep_end
[]
[mc_shear_f]
type = MaterialStdVectorAux
index = 6
property = mc_plastic_yield_function
variable = mc_shear_f
execute_on = timestep_end
[]
[mc_tensile_f]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_yield_function
variable = mc_tensile_f
execute_on = timestep_end
[]
[wp_shear_f]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_yield_function
variable = wp_shear_f
execute_on = timestep_end
[]
[wp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_yield_function
variable = wp_tensile_f
execute_on = timestep_end
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_x
boundary = 'xmin xmax'
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_y
boundary = 'ymin ymax'
value = 0.0
[]
[no_z]
type = DirichletBC
variable = disp_z
boundary = zmin
value = 0.0
[]
[no_wc_x]
type = DirichletBC
variable = wc_x
boundary = 'ymin ymax'
value = 0.0
[]
[no_wc_y]
type = DirichletBC
variable = wc_y
boundary = 'xmin xmax'
value = 0.0
[]
[fix_porepressure]
type = FunctionDirichletBC
variable = porepressure
boundary = 'ymin ymax xmax'
function = ini_pp
[]
[roof_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
pt_vals = '-1E3 1E3'
multipliers = '-1 1'
fluid_phase = 0
flux_function = roof_conductance
boundary = roof
[]
[roof]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = roof
[]
[]
[Functions]
[ini_pp]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0'
symbol_values = '2E3 0.0 1E-5 1E3'
expression = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
[]
[ini_xx]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
[]
[ini_zz]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
[]
[excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '0.5 0 1000.0 1E-9 1 10'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
[]
[density_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '0.5 0 1000.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[]
[roof_conductance]
type = ParsedFunction
symbol_names = 'end_t ymin ymax maxval minval'
symbol_values = '0.5 0 1000.0 1E7 0'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),maxval,minval)'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1 # MPa^-1
[]
[mc_coh_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.99 # MPa
value_residual = 2.01 # MPa
rate = 1.0
[]
[mc_fric]
type = TensorMechanicsHardeningConstant
value = 0.61 # 35deg
[]
[mc_dil]
type = TensorMechanicsHardeningConstant
value = 0.15 # 8deg
[]
[mc_tensile_str_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[]
[mc_compressive_str]
type = TensorMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[]
[wp_coh_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[]
[wp_tan_fric]
type = TensorMechanicsHardeningConstant
value = 0.26 # 15deg
[]
[wp_tan_dil]
type = TensorMechanicsHardeningConstant
value = 0.18 # 10deg
[]
[wp_tensile_str_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[]
[wp_compressive_str_soften]
type = TensorMechanicsHardeningCubic
value_0 = 100
value_residual = 1
internal_limit = 1.0
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E3
density0 = 1000
thermal_expansion = 0
viscosity = 3.5E-17
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity_for_aux]
type = PorousFlowPorosity
at_nodes = false
fluid = true
mechanical = true
ensure_positive = true
porosity_zero = 0.02
solid_bulk = 5.3333E3
[]
[porosity_bulk]
type = PorousFlowPorosity
fluid = true
mechanical = true
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
ensure_positive = true
porosity_zero = 0.02
solid_bulk = 5.3333E3
[]
[porosity_excav]
type = PorousFlowPorosityConst
block = 1
porosity = 1.0
[]
[permeability_bulk]
type = PorousFlowPermeabilityKozenyCarman
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
poroperm_function = kozeny_carman_phi0
k0 = 1E-15
phi0 = 0.02
n = 2
m = 2
[]
[permeability_excav]
type = PorousFlowPermeabilityConst
block = 1
permeability = '0 0 0 0 0 0 0 0 0'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 4
s_res = 0.4
sum_s_res = 0.4
phase = 0
[]
[elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
[]
[elasticity_tensor_1]
type = ComputeLayeredCosseratElasticityTensor
block = 1
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
elasticity_tensor_prefactor = excav_sideways
[]
[strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[]
[ini_stress]
type = ComputeEigenstrainFromInitialStress
eigenstrain_name = ini_stress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
[]
[stress_0]
type = ComputeMultipleInelasticCosseratStress
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
inelastic_models = 'mc wp'
cycle_models = true
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[]
[stress_1]
type = ComputeMultipleInelasticCosseratStress
block = 1
inelastic_models = ''
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[]
[mc]
type = CappedMohrCoulombCosseratStressUpdate
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = mc
tensile_strength = mc_tensile_str_strong_harden
compressive_strength = mc_compressive_str
cohesion = mc_coh_strong_harden
friction_angle = mc_fric
dilation_angle = mc_dil
max_NR_iterations = 100000
smoothing_tol = 0.1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-9 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[]
[wp]
type = CappedWeakPlaneCosseratStressUpdate
warn_about_precision_loss = false
base_name = wp
cohesion = wp_coh_harden
tan_friction_angle = wp_tan_fric
tan_dilation_angle = wp_tan_dil
tensile_strength = wp_tensile_str_harden
compressive_strength = wp_compressive_str_soften
max_NR_iterations = 10000
tip_smoother = 0.05
smoothing_tol = 0.05 # MPa # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0E-3
[]
[undrained_density_0]
type = GenericConstantMaterial
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
prop_names = density
prop_values = 2500
[]
[undrained_density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Postprocessors]
[min_roof_disp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = disp_z
[]
[min_roof_pp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = porepressure
[]
[min_surface_disp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = disp_z
[]
[min_surface_pp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = porepressure
[]
[max_perm_zz]
type = ElementExtremeValue
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
variable = perm_zz
[]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
# best overall
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
# best if you don't have mumps:
#petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' asm 2 lu gmres 200'
# very basic:
#petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' bjacobi gmres 200'
line_search = bt
nl_abs_tol = 1e-3
nl_rel_tol = 1e-5
l_max_its = 200
nl_max_its = 30
start_time = 0.0
dt = 0.0025
end_time = 0.5
[]
[Outputs]
time_step_interval = 1
print_linear_residuals = true
exodus = true
csv = true
console = true
[]
(modules/porous_flow/test/tests/newton_cooling/nc06.i)
# Newton cooling from a bar. 1-phase and heat, steady
[Mesh]
type = GeneratedMesh
dim = 2
nx = 100
ny = 1
xmin = 0
xmax = 100
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pressure temp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Variables]
[pressure]
[]
[temp]
[]
[]
[ICs]
# have to start these reasonably close to their steady-state values
[pressure]
type = FunctionIC
variable = pressure
function = '(2-x/100)*1E6'
[]
[temperature]
type = FunctionIC
variable = temp
function = 100+0.1*x
[]
[]
[Kernels]
[flux]
type = PorousFlowAdvectiveFlux
fluid_component = 0
gravity = '0 0 0'
variable = pressure
[]
[heat_advection]
type = PorousFlowHeatAdvection
gravity = '0 0 0'
variable = temp
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1e6
density0 = 1000
thermal_expansion = 0
viscosity = 1e-3
cv = 1e6
porepressure_coefficient = 0
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
n = 2
phase = 0
[]
[]
[BCs]
[leftp]
type = DirichletBC
variable = pressure
boundary = left
value = 2E6
[]
[leftt]
type = DirichletBC
variable = temp
boundary = left
value = 100
[]
[newtonp]
type = PorousFlowPiecewiseLinearSink
variable = pressure
boundary = right
pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
use_mobility = false
use_relperm = false
fluid_phase = 0
flux_function = 1
[]
[newton]
type = PorousFlowPiecewiseLinearSink
variable = temp
boundary = right
pt_vals = '0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000'
multipliers = '0. 5.6677197748570516e-6 0.000011931518841831313 0.00001885408740732065 0.000026504708864284114 0.000034959953203725676 0.000044304443352900224 0.00005463170211001232 0.00006604508815181467 0.00007865883048198513 0.00009259917167338928 0.00010800563134618119 0.00012503240252705603 0.00014384989486488752 0.00016464644014777016 0.00018763017719085535 0.0002130311349595711 0.00024110353477682344 0.00027212833465544285 0.00030641604122040985 0.00034430981736352295'
use_mobility = false
use_relperm = false
use_internal_energy = true
fluid_phase = 0
flux_function = 1
[]
[]
[VectorPostprocessors]
[porepressure]
type = LineValueSampler
variable = pressure
start_point = '0 0.5 0'
end_point = '100 0.5 0'
sort_by = x
num_points = 11
execute_on = timestep_end
[]
[temperature]
type = LineValueSampler
variable = temp
start_point = '0 0.5 0'
end_point = '100 0.5 0'
sort_by = x
num_points = 11
execute_on = timestep_end
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol '
petsc_options_value = 'gmres asm lu 100 NONZERO 2 1E-8 1E-15'
[]
[]
[Executioner]
type = Steady
solve_type = Newton
[]
[Outputs]
file_base = nc06
execute_on = timestep_end
[along_line]
type = CSV
execute_vector_postprocessors_on = timestep_end
[]
[]
(modules/porous_flow/test/tests/sinks/injection_production_eg.i)
# phase = 0 is liquid phase
# phase = 1 is gas phase
# fluid_component = 0 is water
# fluid_component = 1 is CO2
# Constant rate of CO2 injection into the left boundary
# 1D mesh
# The PorousFlowPiecewiseLinearSinks remove the correct water and CO2 from the right boundary
# Note i take pretty big timesteps here so the system is quite nonlinear
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 20
[]
[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 0 0'
[]
[AuxVariables]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[frac_water_in_liquid]
initial_condition = 1.0
[]
[frac_water_in_gas]
initial_condition = 0.0
[]
[]
[AuxKernels]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = timestep_end
[]
[]
[Variables]
[pwater]
initial_condition = 20E6
[]
[pgas]
initial_condition = 20.1E6
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pwater
[]
[flux0]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pwater
[]
[mass1]
type = PorousFlowMassTimeDerivative
fluid_component = 1
variable = pgas
[]
[flux1]
type = PorousFlowAdvectiveFlux
fluid_component = 1
variable = pgas
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas pwater'
number_fluid_phases = 2
number_fluid_components = 2
[]
[pc]
type = PorousFlowCapillaryPressureVG
alpha = 1E-6
m = 0.6
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[tabulated_water]
type = TabulatedBicubicFluidProperties
fp = true_water
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = water97_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated_11.csv
[]
[true_co2]
type = CO2FluidProperties
[]
[tabulated_co2]
type = TabulatedBicubicFluidProperties
fp = true_co2
temperature_min = 275
pressure_max = 1E8
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = co2_tabulated_11.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = co2_tabulated_11.csv
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = 293.15
[]
[saturation_calculator]
type = PorousFlow2PhasePP
phase0_porepressure = pwater
phase1_porepressure = pgas
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
mass_fraction_vars = 'frac_water_in_liquid frac_water_in_gas'
[]
[water]
type = PorousFlowSingleComponentFluid
fp = tabulated_water
phase = 0
[]
[co2]
type = PorousFlowSingleComponentFluid
fp = tabulated_co2
phase = 1
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
[]
[relperm_water]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
s_res = 0.1
sum_s_res = 0.2
[]
[relperm_gas]
type = PorousFlowRelativePermeabilityBC
nw_phase = true
lambda = 2
s_res = 0.1
sum_s_res = 0.2
phase = 1
[]
[]
[BCs]
[co2_injection]
type = PorousFlowSink
boundary = left
variable = pgas # pgas is associated with the CO2 mass balance (fluid_component = 1 in its Kernels)
flux_function = -1E-2 # negative means a source, rather than a sink
[]
[right_water]
type = PorousFlowPiecewiseLinearSink
boundary = right
# a sink of water, since the Kernels given to pwater are for fluid_component = 0 (the water)
variable = pwater
# this Sink is a function of liquid porepressure
# Also, all the mass_fraction, mobility and relperm are referenced to the liquid phase now
fluid_phase = 0
# Sink strength = (Pwater - 20E6)
pt_vals = '0 1E9'
multipliers = '0 1E9'
PT_shift = 20E6
# multiply Sink strength computed above by mass fraction of water at the boundary
mass_fraction_component = 0
# also multiply Sink strength by mobility of the liquid
use_mobility = true
# also multiply Sink strength by the relperm of the liquid
use_relperm = true
# also multiplly Sink strength by 1/L, where L is the distance to the fixed-porepressure external environment
flux_function = 10 # 1/L
[]
[right_co2]
type = PorousFlowPiecewiseLinearSink
boundary = right
variable = pgas
fluid_phase = 1
pt_vals = '0 1E9'
multipliers = '0 1E9'
PT_shift = 20.1E6
mass_fraction_component = 1
use_mobility = true
use_relperm = true
flux_function = 10 # 1/L
[]
[]
[Preconditioning]
active = 'basic'
[basic]
type = SMP
full = true
petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = 'gmres asm lu NONZERO 2'
[]
[preferred]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
nl_abs_tol = 1E-13
nl_rel_tol = 1E-10
end_time = 1e4
[TimeStepper]
type = IterationAdaptiveDT
dt = 1E4
growth_factor = 1.1
[]
[]
[VectorPostprocessors]
[pps]
type = LineValueSampler
warn_discontinuous_face_values = false
start_point = '0 0 0'
end_point = '20 0 0'
num_points = 20
sort_by = x
variable = 'pgas pwater saturation_gas'
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/numerical_diffusion/fully_saturated_action.i)
# Using the fully-saturated action, which does mass lumping but no upwinding
[Mesh]
type = GeneratedMesh
dim = 1
nx = 100
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
[]
[tracer]
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '1 - x'
[]
[tracer]
type = FunctionIC
variable = tracer
function = 'if(x<0.1,0,if(x>0.3,0,1))'
[]
[]
[PorousFlowFullySaturated]
porepressure = porepressure
coupling_type = Hydro
gravity = '0 0 0'
fp = the_simple_fluid
mass_fraction_vars = tracer
stabilization = none
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1
boundary = left
[]
[no_tracer_on_left]
type = DirichletBC
variable = tracer
value = 0
boundary = left
[]
[remove_component_1]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 1
use_mobility = true
flux_function = 1E3
[]
[remove_component_0]
type = PorousFlowPiecewiseLinearSink
variable = tracer
boundary = right
fluid_phase = 0
pt_vals = '0 1E3'
multipliers = '0 1E3'
mass_fraction_component = 0
use_mobility = true
flux_function = 1E3
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
thermal_expansion = 0
viscosity = 1.0
density0 = 1000.0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-2 0 0 0 1E-2 0 0 0 1E-2'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[VectorPostprocessors]
[tracer]
type = LineValueSampler
start_point = '0 0 0'
end_point = '1 0 0'
num_points = 101
sort_by = x
variable = tracer
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 6
dt = 6E-1
nl_abs_tol = 1E-8
timestep_tolerance = 1E-3
[]
[Outputs]
[out]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/test/tests/newton_cooling/nc08.i)
# Newton cooling from a bar. 1-phase ideal fluid and heat, steady
[Mesh]
type = GeneratedMesh
dim = 2
nx = 100
ny = 1
xmin = 0
xmax = 100
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pressure temp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Variables]
[pressure]
[]
[temp]
[]
[]
[ICs]
# have to start these reasonably close to their steady-state values
[pressure]
type = FunctionIC
variable = pressure
function = '200-0.5*x'
[]
[temperature]
type = FunctionIC
variable = temp
function = 180+0.1*x
[]
[]
[Kernels]
[flux]
type = PorousFlowAdvectiveFlux
fluid_component = 0
gravity = '0 0 0'
variable = pressure
[]
[heat_advection]
type = PorousFlowHeatAdvection
gravity = '0 0 0'
variable = temp
[]
[]
[FluidProperties]
[idealgas]
type = IdealGasFluidProperties
molar_mass = 1.4
gamma = 1.2
mu = 1.2
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[dens0]
type = PorousFlowSingleComponentFluid
fp = idealgas
phase = 0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.1 0 0 0 1.1 0 0 0 1.1'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey # irrelevant in this fully-saturated situation
n = 2
phase = 0
[]
[]
[BCs]
[leftp]
type = DirichletBC
variable = pressure
boundary = left
value = 200
[]
[leftt]
type = DirichletBC
variable = temp
boundary = left
value = 180
[]
[newtonp]
type = PorousFlowPiecewiseLinearSink
variable = pressure
boundary = right
pt_vals = '-200 0 200'
multipliers = '-200 0 200'
use_mobility = true
use_relperm = true
fluid_phase = 0
flux_function = 0.005 # 1/2/L
[]
[newtont]
type = PorousFlowPiecewiseLinearSink
variable = temp
boundary = right
pt_vals = '-200 0 200'
multipliers = '-200 0 200'
use_mobility = true
use_relperm = true
use_enthalpy = true
fluid_phase = 0
flux_function = 0.005 # 1/2/L
[]
[]
[VectorPostprocessors]
[porepressure]
type = LineValueSampler
variable = pressure
start_point = '0 0.5 0'
end_point = '100 0.5 0'
sort_by = x
num_points = 11
execute_on = timestep_end
[]
[temperature]
type = LineValueSampler
variable = temp
start_point = '0 0.5 0'
end_point = '100 0.5 0'
sort_by = x
num_points = 11
execute_on = timestep_end
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
[]
[]
[Executioner]
type = Steady
solve_type = Newton
nl_rel_tol = 1E-10
nl_abs_tol = 1E-15
[]
[Outputs]
file_base = nc08
execute_on = timestep_end
[along_line]
type = CSV
execute_vector_postprocessors_on = timestep_end
[]
[]