- postprocessorThe name of the postprocessor that computes the dt
C++ Type:PostprocessorName
Description:The name of the postprocessor that computes the dt
PostprocessorDT
under construction:Undocumented Class
The PostprocessorDT has not been documented. The content listed below should be used as a starting point for documenting the class, which includes the typical automatic documentation associated with a MooseObject; however, what is contained is ultimately determined by what is necessary to make the documentation clear for users.
# PostprocessorDT
!syntax description /Executioner/TimeStepper/PostprocessorDT
## Overview
!! Replace these lines with information regarding the PostprocessorDT object.
## Example Input File Syntax
!! Describe and include an example of how to use the PostprocessorDT object.
!syntax parameters /Executioner/TimeStepper/PostprocessorDT
!syntax inputs /Executioner/TimeStepper/PostprocessorDT
!syntax children /Executioner/TimeStepper/PostprocessorDT
!syntax description /Executioner/TimeStepper/PostprocessorDT
Input Parameters
- cutback_factor_at_failure0.5Factor to apply to timestep if it a time step fails to convergence.
Default:0.5
C++ Type:double
Options:
Description:Factor to apply to timestep if it a time step fails to convergence.
- dtInitial value of dt
C++ Type:double
Options:
Description:Initial value of dt
- factor0Add a factor to the supplied postprocessor value.
Default:0
C++ Type:double
Options:
Description:Add a factor to the supplied postprocessor value.
- reset_dtFalseUse when restarting a calculation to force a change in dt.
Default:False
C++ Type:bool
Options:
Description:Use when restarting a calculation to force a change in dt.
- scale1Multiple scale and supplied postprocessor value.
Default:1
C++ Type:double
Options:
Description:Multiple scale and supplied postprocessor value.
Optional Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector
Options:
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Options:
Description:Set the enabled status of the MooseObject.
Advanced Parameters
Input Files
- modules/level_set/test/tests/reinitialization/master.i
- modules/combined/test/tests/poro_mechanics/borehole_highres.i
- modules/level_set/examples/vortex/vortex_supg.i
- modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated.i
- modules/porous_flow/test/tests/poro_elasticity/terzaghi.i
- modules/level_set/examples/rotating_circle/circle_rotate.i
- modules/porous_flow/test/tests/poro_elasticity/terzaghi_fully_saturated_volume.i
- modules/porous_flow/test/tests/poro_elasticity/terzaghi_basicthm.i
- modules/porous_flow/test/tests/poro_elasticity/terzaghi_constM.i
- modules/combined/test/tests/poro_mechanics/borehole_lowres.i
- modules/porous_flow/test/tests/poro_elasticity/mandel_constM.i
- modules/level_set/examples/rotating_circle/circle_rotate_supg.i
- modules/level_set/examples/vortex/vortex_reinit.i
- modules/combined/test/tests/poro_mechanics/terzaghi.i
- modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated_volume.i
- modules/level_set/examples/circle/circle_16.i
- modules/level_set/examples/rotating_circle/circle_rotate_master.i
- modules/porous_flow/test/tests/poro_elasticity/mandel_basicthm.i
- test/tests/time_steppers/postprocessor_dt/postprocessor_dt.i
- modules/porous_flow/test/tests/poro_elasticity/mandel.i
- modules/combined/test/tests/poro_mechanics/mandel.i
- modules/level_set/examples/vortex/vortex.i
modules/level_set/test/tests/reinitialization/master.i
[Mesh]
type = GeneratedMesh
dim = 2
xmin = 0
xmax = 1
ymin = 0
ymax = 1
nx = 8
ny = 8
uniform_refine = 3 #1/64
[]
[AuxVariables]
[./vel_x]
initial_condition = 1
[../]
[./vel_y]
initial_condition = 1
[../]
[]
[Variables]
[./phi]
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.05
center = '0.5 0.5 0'
radius = 0.15
[../]
[]
[BCs]
[./Periodic]
[./all]
variable = phi
auto_direction = 'x y'
[../]
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[../]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 1
nl_rel_tol = 1e-12
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 1
[../]
[]
[MultiApps]
[./reinit]
type = LevelSetReinitializationMultiApp
input_files = 'reinit.i'
execute_on = 'timestep_end'
[../]
[]
[Transfers]
[./to_sub]
type = MultiAppCopyTransfer
variable = phi
source_variable = phi
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[./to_sub_init]
type = MultiAppCopyTransfer
variable = phi_0
source_variable = phi
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[./from_sub]
type = MultiAppCopyTransfer
variable = phi
source_variable = phi
direction = from_multiapp
multi_app = reinit
execute_on = timestep_end
[../]
[]
[Outputs]
exodus = true
csv = true
[]
modules/combined/test/tests/poro_mechanics/borehole_highres.i
# Poroelastic response of a borehole.
#
# HIGHRES VERSION: this version gives good agreement with the analytical solution, but it takes a while so is a "heavy" test
#
# A fully-saturated medium contains a fluid with a homogeneous porepressure,
# but an anisitropic insitu stress. A infinitely-long borehole aligned with
# the $$z$$ axis is instanteously excavated. The borehole boundary is
# stress-free and allowed to freely drain. This problem is analysed using
# plane-strain conditions (no $$z$$ displacement).
#
# The solution in Laplace space is found in E Detournay and AHD Cheng "Poroelastic response of a borehole in a non-hydrostatic stress field". International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 25 (1988) 171-182. In the small-time limit, the Laplace transforms may be performed. There is one typo in the paper. Equation (A4)'s final term should be -(a/r)\sqrt(4ct/(a^2\pi)), and not +(a/r)\sqrt(4ct/(a^2\pi)).
#
# Because realistic parameters are chosen (below),
# the residual for porepressure is much smaller than
# the residuals for the displacements. Therefore the
# scaling parameter is chosen. Also note that the
# insitu stresses are effective stresses, not total
# stresses, but the solution in the above paper is
# expressed in terms of total stresses.
#
# Here are the problem's parameters, and their values:
# Borehole radius. a = 1
# Rock's Lame lambda. la = 0.5E9
# Rock's Lame mu, which is also the Rock's shear modulus. mu = G = 1.5E9
# Rock bulk modulus. K = la + 2*mu/3 = 1.5E9
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.125
# Rock bulk compliance. 1/K = 0.66666666E-9
# Fluid bulk modulus. Kf = 0.7171315E9
# Fluid bulk compliance. 1/Kf = 1.39444444E-9
# Rock initial porosity. phi0 = 0.3
# Biot coefficient. alpha = 0.65
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 2E9
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.345E9
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.2364
# Skempton coefficient. B = alpha*M/Ku = 0.554
# Fluid mobility (rock permeability/fluid viscosity). k = 1E-12
[Mesh]
type = FileMesh
file = borehole_highres_input.e
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
block = 1
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
scaling = 1E9 # Notice the scaling, to make porepressure's kernels roughly of same magnitude as disp's kernels
[../]
[]
[GlobalParams]
volumetric_locking_correction=true
[]
[ICs]
[./initial_p]
type = ConstantIC
variable = porepressure
value = 1E6
[../]
[]
[BCs]
[./fixed_outer_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = outer
[../]
[./fixed_outer_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = outer
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'zmin zmax'
[../]
[./borehole_wall]
type = DirichletBC
variable = porepressure
value = 0
boundary = bh_wall
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./tot_yy]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./tot_yy]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_yy
function = 'stress_yy-0.65*porepressure'
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PoroMechanicsCoupling
variable = disp_x
component = 0
[../]
[./poro_y]
type = PoroMechanicsCoupling
variable = disp_y
component = 1
[../]
[./poro_z]
type = PoroMechanicsCoupling
variable = disp_z
component = 2
[../]
[./poro_timederiv]
type = PoroFullSatTimeDerivative
variable = porepressure
[../]
[./darcy_flow]
type = CoefDiffusion
variable = porepressure
coef = 1E-12
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5E9 1.5E9'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*1.5/3 = 1.5E9
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeFiniteStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '-1.35E6 0 0 0 -3.35E6 0 0 0 0' # remember this is the effective stress
eigenstrain_name = ini_stress
[../]
[./no_plasticity]
type = ComputeFiniteStrainElasticStress
[../]
[./poro_material]
type = PoroFullSatMaterial
porosity0 = 0.3
biot_coefficient = 0.65
solid_bulk_compliance = 0.6666666666667E-9
fluid_bulk_compliance = 1.3944444444444E-9
constant_porosity = false
[../]
[]
[Postprocessors]
[./p00]
type = PointValue
variable = porepressure
point = '1.00 0 0'
outputs = csv_p
[../]
[./p01]
type = PointValue
variable = porepressure
point = '1.01 0 0'
outputs = csv_p
[../]
[./p02]
type = PointValue
variable = porepressure
point = '1.02 0 0'
outputs = csv_p
[../]
[./p03]
type = PointValue
variable = porepressure
point = '1.03 0 0'
outputs = csv_p
[../]
[./p04]
type = PointValue
variable = porepressure
point = '1.04 0 0'
outputs = csv_p
[../]
[./p05]
type = PointValue
variable = porepressure
point = '1.05 0 0'
outputs = csv_p
[../]
[./p06]
type = PointValue
variable = porepressure
point = '1.06 0 0'
outputs = csv_p
[../]
[./p07]
type = PointValue
variable = porepressure
point = '1.07 0 0'
outputs = csv_p
[../]
[./p08]
type = PointValue
variable = porepressure
point = '1.08 0 0'
outputs = csv_p
[../]
[./p09]
type = PointValue
variable = porepressure
point = '1.09 0 0'
outputs = csv_p
[../]
[./p10]
type = PointValue
variable = porepressure
point = '1.10 0 0'
outputs = csv_p
[../]
[./p11]
type = PointValue
variable = porepressure
point = '1.11 0 0'
outputs = csv_p
[../]
[./p12]
type = PointValue
variable = porepressure
point = '1.12 0 0'
outputs = csv_p
[../]
[./p13]
type = PointValue
variable = porepressure
point = '1.13 0 0'
outputs = csv_p
[../]
[./p14]
type = PointValue
variable = porepressure
point = '1.14 0 0'
outputs = csv_p
[../]
[./p15]
type = PointValue
variable = porepressure
point = '1.15 0 0'
outputs = csv_p
[../]
[./p16]
type = PointValue
variable = porepressure
point = '1.16 0 0'
outputs = csv_p
[../]
[./p17]
type = PointValue
variable = porepressure
point = '1.17 0 0'
outputs = csv_p
[../]
[./p18]
type = PointValue
variable = porepressure
point = '1.18 0 0'
outputs = csv_p
[../]
[./p19]
type = PointValue
variable = porepressure
point = '1.19 0 0'
outputs = csv_p
[../]
[./p20]
type = PointValue
variable = porepressure
point = '1.20 0 0'
outputs = csv_p
[../]
[./p21]
type = PointValue
variable = porepressure
point = '1.21 0 0'
outputs = csv_p
[../]
[./p22]
type = PointValue
variable = porepressure
point = '1.22 0 0'
outputs = csv_p
[../]
[./p23]
type = PointValue
variable = porepressure
point = '1.23 0 0'
outputs = csv_p
[../]
[./p24]
type = PointValue
variable = porepressure
point = '1.24 0 0'
outputs = csv_p
[../]
[./p25]
type = PointValue
variable = porepressure
point = '1.25 0 0'
outputs = csv_p
[../]
[./s00]
type = PointValue
variable = disp_x
point = '1.00 0 0'
outputs = csv_s
[../]
[./s01]
type = PointValue
variable = disp_x
point = '1.01 0 0'
outputs = csv_s
[../]
[./s02]
type = PointValue
variable = disp_x
point = '1.02 0 0'
outputs = csv_s
[../]
[./s03]
type = PointValue
variable = disp_x
point = '1.03 0 0'
outputs = csv_s
[../]
[./s04]
type = PointValue
variable = disp_x
point = '1.04 0 0'
outputs = csv_s
[../]
[./s05]
type = PointValue
variable = disp_x
point = '1.05 0 0'
outputs = csv_s
[../]
[./s06]
type = PointValue
variable = disp_x
point = '1.06 0 0'
outputs = csv_s
[../]
[./s07]
type = PointValue
variable = disp_x
point = '1.07 0 0'
outputs = csv_s
[../]
[./s08]
type = PointValue
variable = disp_x
point = '1.08 0 0'
outputs = csv_s
[../]
[./s09]
type = PointValue
variable = disp_x
point = '1.09 0 0'
outputs = csv_s
[../]
[./s10]
type = PointValue
variable = disp_x
point = '1.10 0 0'
outputs = csv_s
[../]
[./s11]
type = PointValue
variable = disp_x
point = '1.11 0 0'
outputs = csv_s
[../]
[./s12]
type = PointValue
variable = disp_x
point = '1.12 0 0'
outputs = csv_s
[../]
[./s13]
type = PointValue
variable = disp_x
point = '1.13 0 0'
outputs = csv_s
[../]
[./s14]
type = PointValue
variable = disp_x
point = '1.14 0 0'
outputs = csv_s
[../]
[./s15]
type = PointValue
variable = disp_x
point = '1.15 0 0'
outputs = csv_s
[../]
[./s16]
type = PointValue
variable = disp_x
point = '1.16 0 0'
outputs = csv_s
[../]
[./s17]
type = PointValue
variable = disp_x
point = '1.17 0 0'
outputs = csv_s
[../]
[./s18]
type = PointValue
variable = disp_x
point = '1.18 0 0'
outputs = csv_s
[../]
[./s19]
type = PointValue
variable = disp_x
point = '1.19 0 0'
outputs = csv_s
[../]
[./s20]
type = PointValue
variable = disp_x
point = '1.20 0 0'
outputs = csv_s
[../]
[./s21]
type = PointValue
variable = disp_x
point = '1.21 0 0'
outputs = csv_s
[../]
[./s22]
type = PointValue
variable = disp_x
point = '1.22 0 0'
outputs = csv_s
[../]
[./s23]
type = PointValue
variable = disp_x
point = '1.23 0 0'
outputs = csv_s
[../]
[./s24]
type = PointValue
variable = disp_x
point = '1.24 0 0'
outputs = csv_s
[../]
[./s25]
type = PointValue
variable = disp_x
point = '1.25 0 0'
outputs = csv_s
[../]
[./t00]
type = PointValue
variable = tot_yy
point = '1.00 0 0'
outputs = csv_t
[../]
[./t01]
type = PointValue
variable = tot_yy
point = '1.01 0 0'
outputs = csv_t
[../]
[./t02]
type = PointValue
variable = tot_yy
point = '1.02 0 0'
outputs = csv_t
[../]
[./t03]
type = PointValue
variable = tot_yy
point = '1.03 0 0'
outputs = csv_t
[../]
[./t04]
type = PointValue
variable = tot_yy
point = '1.04 0 0'
outputs = csv_t
[../]
[./t05]
type = PointValue
variable = tot_yy
point = '1.05 0 0'
outputs = csv_t
[../]
[./t06]
type = PointValue
variable = tot_yy
point = '1.06 0 0'
outputs = csv_t
[../]
[./t07]
type = PointValue
variable = tot_yy
point = '1.07 0 0'
outputs = csv_t
[../]
[./t08]
type = PointValue
variable = tot_yy
point = '1.08 0 0'
outputs = csv_t
[../]
[./t09]
type = PointValue
variable = tot_yy
point = '1.09 0 0'
outputs = csv_t
[../]
[./t10]
type = PointValue
variable = tot_yy
point = '1.10 0 0'
outputs = csv_t
[../]
[./t11]
type = PointValue
variable = tot_yy
point = '1.11 0 0'
outputs = csv_t
[../]
[./t12]
type = PointValue
variable = tot_yy
point = '1.12 0 0'
outputs = csv_t
[../]
[./t13]
type = PointValue
variable = tot_yy
point = '1.13 0 0'
outputs = csv_t
[../]
[./t14]
type = PointValue
variable = tot_yy
point = '1.14 0 0'
outputs = csv_t
[../]
[./t15]
type = PointValue
variable = tot_yy
point = '1.15 0 0'
outputs = csv_t
[../]
[./t16]
type = PointValue
variable = tot_yy
point = '1.16 0 0'
outputs = csv_t
[../]
[./t17]
type = PointValue
variable = tot_yy
point = '1.17 0 0'
outputs = csv_t
[../]
[./t18]
type = PointValue
variable = tot_yy
point = '1.18 0 0'
outputs = csv_t
[../]
[./t19]
type = PointValue
variable = tot_yy
point = '1.19 0 0'
outputs = csv_t
[../]
[./t20]
type = PointValue
variable = tot_yy
point = '1.20 0 0'
outputs = csv_t
[../]
[./t21]
type = PointValue
variable = tot_yy
point = '1.21 0 0'
outputs = csv_t
[../]
[./t22]
type = PointValue
variable = tot_yy
point = '1.22 0 0'
outputs = csv_t
[../]
[./t23]
type = PointValue
variable = tot_yy
point = '1.23 0 0'
outputs = csv_t
[../]
[./t24]
type = PointValue
variable = tot_yy
point = '1.24 0 0'
outputs = csv_t
[../]
[./t25]
type = PointValue
variable = tot_yy
point = '1.25 0 0'
outputs = csv_t
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = 2*t
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options = '-snes_monitor -snes_linesearch_monitor'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm 1E0 1E-10 200 500 lu NONZERO'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.3
dt = 0.1
#[./TimeStepper]
# type = PostprocessorDT
# postprocessor = dt
# dt = 0.003
#[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = borehole_highres
exodus = true
sync_times = '0.003 0.3'
[./csv_p]
file_base = borehole_highres_p
type = CSV
[../]
[./csv_s]
file_base = borehole_highres_s
type = CSV
[../]
[./csv_t]
file_base = borehole_highres_t
type = CSV
[../]
[]
modules/level_set/examples/vortex/vortex_supg.i
[Mesh]
type = GeneratedMesh
dim = 2
xmax = 1
ymax = 1
nx = 16
ny = 16
uniform_refine = 2
elem_type = QUAD9
[]
[AuxVariables]
[./vel_x]
family = LAGRANGE
order = FIRST
[../]
[./vel_y]
family = LAGRANGE
order = FIRST
[../]
[]
[AuxKernels]
[./vel_x]
type = FunctionAux
function = vel_x
variable = vel_x
execute_on = 'initial timestep_begin'
[../]
[./vel_y]
type = FunctionAux
function = vel_y
variable = vel_y
execute_on = 'initial timestep_begin'
[../]
[]
[Variables]
[./phi]
family = LAGRANGE
order = FIRST
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.01184
center = '0.5 0.75 0'
radius = 0.15
[../]
[./vel_x]
type = LevelSetOlssonVortex
component = x
reverse_time = 2
[../]
[./vel_y]
type = LevelSetOlssonVortex
component = y
reverse_time = 2
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[./advection_supg]
type = LevelSetAdvectionSUPG
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[./time_supg]
type = LevelSetTimeDerivativeSUPG
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[../]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial timestep_end'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 2
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[../]
[]
[Outputs]
csv = true
exodus = true
[]
modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated.i
# Mandel's problem of consolodation of a drained medium
# Using the FullySaturatedDarcyBase and FullySaturatedMassTimeDerivative kernels
# Note the use of consistent_with_displaced_mesh = false in the calculation of volumetric strain
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[../]
[./roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[../]
[./xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[../]
[./top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[../]
[]
[Functions]
[./top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./tot_force]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./tot_force]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
function = '-stress_yy+0.6*porepressure'
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[../]
[./mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
biot_coefficient = 0.6
coupling_type = HydroMechanical
variable = porepressure
[../]
[./flux]
type = PorousFlowFullySaturatedDarcyBase
variable = porepressure
gravity = '0 0 0'
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 1
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./eff_fluid_pressure_qp]
type = PorousFlowEffectiveFluidPressure
[../]
[./vol_strain]
type = PorousFlowVolumetricStrain
consistent_with_displaced_mesh = false
[../]
[./ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = porepressure
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid_qp]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[../]
[./porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[../]
[./biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
solid_bulk_compliance = 1
fluid_bulk_modulus = 8
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[../]
[./p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[../]
[./p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[../]
[./p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[../]
[./p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[../]
[./p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[../]
[./p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[../]
[./p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[../]
[./p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[../]
[./p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[../]
[./p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[../]
[./xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[../]
[./ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[../]
[./total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel_fully_saturated
[./csv]
interval = 3
type = CSV
[../]
[]
modules/porous_flow/test/tests/poro_elasticity/terzaghi.i
# Terzaghi's problem of consolodation of a drained medium
#
# A saturated soil sample sits in a bath of water.
# It is constrained on its sides, and bottom.
# Its sides and bottom are also impermeable.
# Initially it is unstressed.
# A normal stress, q, is applied to the soil's top.
# The soil then slowly compresses as water is squeezed
# out from the sample from its top (the top BC for
# the porepressure is porepressure = 0).
#
# See, for example. Section 2.2 of the online manuscript
# Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013
# but note that the "sigma" in that paper is the negative
# of the stress in TensorMechanics
#
# Here are the problem's parameters, and their values:
# Soil height. h = 10
# Soil's Lame lambda. la = 2
# Soil's Lame mu, which is also the Soil's shear modulus. mu = 3
# Soil bulk modulus. K = la + 2*mu/3 = 4
# Soil confined compressibility. m = 1/(K + 4mu/3) = 0.125
# Soil bulk compliance. 1/K = 0.25
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Soil initial storativity, which is the reciprocal of the initial Biot modulus. S = phi0/Kf + (alpha - phi0)(1 - alpha)/K = 0.0625
# Consolidation coefficient. c = k/(S + alpha^2 m) = 13.95348837
# Normal stress on top. q = 1
# Initial porepressure, resulting from instantaneous application of q, assuming corresponding instantaneous increase of porepressure (Note that this is calculated by MOOSE: we only need it for the analytical solution). p0 = alpha*m*q/(S + alpha^2 m) = 0.69767442
# Initial vertical displacement (down is positive), resulting from instantaneous application of q (Note this is calculated by MOOSE: we only need it for the analytical solution). uz0 = q*m*h*S/(S + alpha^2 m)
# Final vertical displacement (down in positive) (Note this is calculated by MOOSE: we only need it for the analytical solution). uzinf = q*m*h
#
# The solution for porepressure is
# P = 4*p0/\pi \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{2k-1} \cos ((2k-1)\pi z/(2h)) \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
# This series converges very slowly for ct/h^2 small, so in that domain
# P = p0 erf( (1-(z/h))/(2 \sqrt(ct/h^2)) )
#
# The degree of consolidation is defined as
# U = (uz - uz0)/(uzinf - uz0)
# where uz0 and uzinf are defined above, and
# uz = the vertical displacement of the top (down is positive)
# U = 1 - (8/\pi^2)\sum_{k=1}^{\infty} \frac{1}{(2k-1)^2} \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
#
# FINAL NOTE: The above solution assumes constant Biot Modulus.
# In porous_flow this is not true. Therefore the solution is
# a little different than in the paper. This test was therefore
# validated against MOOSE's poromechanics, which can choose either
# a constant Biot Modulus (which has been shown to agree with
# the analytic solution), or a non-constant Biot Modulus (which
# gives the same results as porous_flow).
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 10
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = 0
zmax = 10
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[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.8
alpha = 1
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[../]
[./confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[../]
[./basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[../]
[./topdrained]
type = DirichletBC
variable = porepressure
value = 0
boundary = front
[../]
[./topload]
type = NeumannBC
variable = disp_z
value = -1
boundary = front
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[../]
[./poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[../]
[./mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[../]
[./flux]
type = PorousFlowAdvectiveFlux
variable = porepressure
gravity = '0 0 0'
fluid_component = 0
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 0.96
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '2 3'
# bulk modulus is lambda + 2*mu/3 = 2 + 2*3/3 = 4
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./eff_fluid_pressure_qp]
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]
type = PorousFlowPorosity
fluid = true
mechanical = true
ensure_positive = false
porosity_zero = 0.1
biot_coefficient = 0.6
solid_bulk = 4
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0 # unimportant in this fully-saturated situation
phase = 0
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
use_displaced_mesh = false
[../]
[./p1]
type = PointValue
outputs = csv
point = '0 0 1'
variable = porepressure
use_displaced_mesh = false
[../]
[./p2]
type = PointValue
outputs = csv
point = '0 0 2'
variable = porepressure
use_displaced_mesh = false
[../]
[./p3]
type = PointValue
outputs = csv
point = '0 0 3'
variable = porepressure
use_displaced_mesh = false
[../]
[./p4]
type = PointValue
outputs = csv
point = '0 0 4'
variable = porepressure
use_displaced_mesh = false
[../]
[./p5]
type = PointValue
outputs = csv
point = '0 0 5'
variable = porepressure
use_displaced_mesh = false
[../]
[./p6]
type = PointValue
outputs = csv
point = '0 0 6'
variable = porepressure
use_displaced_mesh = false
[../]
[./p7]
type = PointValue
outputs = csv
point = '0 0 7'
variable = porepressure
use_displaced_mesh = false
[../]
[./p8]
type = PointValue
outputs = csv
point = '0 0 8'
variable = porepressure
use_displaced_mesh = false
[../]
[./p9]
type = PointValue
outputs = csv
point = '0 0 9'
variable = porepressure
use_displaced_mesh = false
[../]
[./p99]
type = PointValue
outputs = csv
point = '0 0 10'
variable = porepressure
use_displaced_mesh = false
[../]
[./zdisp]
type = PointValue
outputs = csv
point = '0 0 10'
variable = disp_z
use_displaced_mesh = false
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.5*t<0.1,0.5*t,0.1)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.0001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = terzaghi
[./csv]
type = CSV
[../]
[]
modules/level_set/examples/rotating_circle/circle_rotate.i
[Mesh]
type = GeneratedMesh
dim = 2
xmin = -1
xmax = 1
ymin = -1
ymax = 1
nx = 32
ny = 32
uniform_refine = 2
[]
[AuxVariables]
[./vel_x]
[../]
[./vel_y]
[../]
[]
[AuxKernels]
[./vel_x]
type = FunctionAux
function = 4*y
variable = vel_x
execute_on = initial
[../]
[./vel_y]
type = FunctionAux
function = -4*x
variable = vel_y
execute_on = initial
[../]
[]
[Variables]
[./phi]
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.03
center = '0 0.5 0'
radius = 0.15
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[../]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial' #timestep_end'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 1.570796
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[../]
[]
[Outputs]
csv = true
exodus = true
[]
modules/porous_flow/test/tests/poro_elasticity/terzaghi_fully_saturated_volume.i
# Terzaghi's problem of consolodation of a drained medium
# The FullySaturated Kernels are used, with multiply_by_density = false
# so that this becomes a linear problem with constant Biot Modulus
# Also, since the FullySaturated Kernels are used, we have to
# use consistent_with_displaced_mesh = false in the calculation of volumetric strain
#
# A saturated soil sample sits in a bath of water.
# It is constrained on its sides, and bottom.
# Its sides and bottom are also impermeable.
# Initially it is unstressed.
# A normal stress, q, is applied to the soil's top.
# The soil then slowly compresses as water is squeezed
# out from the sample from its top (the top BC for
# the porepressure is porepressure = 0).
#
# See, for example. Section 2.2 of the online manuscript
# Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013
# but note that the "sigma" in that paper is the negative
# of the stress in TensorMechanics
#
# Here are the problem's parameters, and their values:
# Soil height. h = 10
# Soil's Lame lambda. la = 2
# Soil's Lame mu, which is also the Soil's shear modulus. mu = 3
# Soil bulk modulus. K = la + 2*mu/3 = 4
# Soil confined compressibility. m = 1/(K + 4mu/3) = 0.125
# Soil bulk compliance. 1/K = 0.25
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Soil initial storativity, which is the reciprocal of the initial Biot modulus. S = phi0/Kf + (alpha - phi0)(1 - alpha)/K = 0.0625
# Consolidation coefficient. c = k/(S + alpha^2 m) = 13.95348837
# Normal stress on top. q = 1
# Initial porepressure, resulting from instantaneous application of q, assuming corresponding instantaneous increase of porepressure (Note that this is calculated by MOOSE: we only need it for the analytical solution). p0 = alpha*m*q/(S + alpha^2 m) = 0.69767442
# Initial vertical displacement (down is positive), resulting from instantaneous application of q (Note this is calculated by MOOSE: we only need it for the analytical solution). uz0 = q*m*h*S/(S + alpha^2 m)
# Final vertical displacement (down in positive) (Note this is calculated by MOOSE: we only need it for the analytical solution). uzinf = q*m*h
#
# The solution for porepressure is
# P = 4*p0/\pi \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{2k-1} \cos ((2k-1)\pi z/(2h)) \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
# This series converges very slowly for ct/h^2 small, so in that domain
# P = p0 erf( (1-(z/h))/(2 \sqrt(ct/h^2)) )
#
# The degree of consolidation is defined as
# U = (uz - uz0)/(uzinf - uz0)
# where uz0 and uzinf are defined above, and
# uz = the vertical displacement of the top (down is positive)
# U = 1 - (8/\pi^2)\sum_{k=1}^{\infty} \frac{1}{(2k-1)^2} \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 10
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = 0
zmax = 10
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[../]
[./confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[../]
[./basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[../]
[./topdrained]
type = DirichletBC
variable = porepressure
value = 0
boundary = front
[../]
[./topload]
type = NeumannBC
variable = disp_z
value = -1
boundary = front
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[../]
[./mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
coupling_type = HydroMechanical
biot_coefficient = 0.6
multiply_by_density = false
variable = porepressure
[../]
[./flux]
type = PorousFlowFullySaturatedDarcyBase
multiply_by_density = false
variable = porepressure
gravity = '0 0 0'
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 0.96
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '2 3'
# bulk modulus is lambda + 2*mu/3 = 2 + 2*3/3 = 4
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./eff_fluid_pressure_qp]
type = PorousFlowEffectiveFluidPressure
[../]
[./vol_strain]
type = PorousFlowVolumetricStrain
consistent_with_displaced_mesh = false
[../]
[./ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = porepressure
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid_qp]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[../]
[./porosity]
type = PorousFlowPorosityConst # only the initial value of this is used
porosity = 0.1
[../]
[./biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
fluid_bulk_modulus = 8
solid_bulk_compliance = 0.25
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
use_displaced_mesh = false
[../]
[./p1]
type = PointValue
outputs = csv
point = '0 0 1'
variable = porepressure
use_displaced_mesh = false
[../]
[./p2]
type = PointValue
outputs = csv
point = '0 0 2'
variable = porepressure
use_displaced_mesh = false
[../]
[./p3]
type = PointValue
outputs = csv
point = '0 0 3'
variable = porepressure
use_displaced_mesh = false
[../]
[./p4]
type = PointValue
outputs = csv
point = '0 0 4'
variable = porepressure
use_displaced_mesh = false
[../]
[./p5]
type = PointValue
outputs = csv
point = '0 0 5'
variable = porepressure
use_displaced_mesh = false
[../]
[./p6]
type = PointValue
outputs = csv
point = '0 0 6'
variable = porepressure
use_displaced_mesh = false
[../]
[./p7]
type = PointValue
outputs = csv
point = '0 0 7'
variable = porepressure
use_displaced_mesh = false
[../]
[./p8]
type = PointValue
outputs = csv
point = '0 0 8'
variable = porepressure
use_displaced_mesh = false
[../]
[./p9]
type = PointValue
outputs = csv
point = '0 0 9'
variable = porepressure
use_displaced_mesh = false
[../]
[./p99]
type = PointValue
outputs = csv
point = '0 0 10'
variable = porepressure
use_displaced_mesh = false
[../]
[./zdisp]
type = PointValue
outputs = csv
point = '0 0 10'
variable = disp_z
use_displaced_mesh = false
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.5*t<0.1,0.5*t,0.1)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.0001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = terzaghi_fully_saturated_volume
[./csv]
type = CSV
[../]
[]
modules/porous_flow/test/tests/poro_elasticity/terzaghi_basicthm.i
# Using a BasicTHM action
# Terzaghi's problem of consolodation of a drained medium
# The FullySaturated Kernels are used, with multiply_by_density = false
# so that this becomes a linear problem with constant Biot Modulus
# Also, since the FullySaturated Kernels are used, we have to
# use consistent_with_displaced_mesh = false in the calculation of volumetric strain
#
# A saturated soil sample sits in a bath of water.
# It is constrained on its sides, and bottom.
# Its sides and bottom are also impermeable.
# Initially it is unstressed.
# A normal stress, q, is applied to the soil's top.
# The soil then slowly compresses as water is squeezed
# out from the sample from its top (the top BC for
# the porepressure is porepressure = 0).
#
# See, for example. Section 2.2 of the online manuscript
# Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013
# but note that the "sigma" in that paper is the negative
# of the stress in TensorMechanics
#
# Here are the problem's parameters, and their values:
# Soil height. h = 10
# Soil's Lame lambda. la = 2
# Soil's Lame mu, which is also the Soil's shear modulus. mu = 3
# Soil bulk modulus. K = la + 2*mu/3 = 4
# Soil confined compressibility. m = 1/(K + 4mu/3) = 0.125
# Soil bulk compliance. 1/K = 0.25
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Soil initial storativity, which is the reciprocal of the initial Biot modulus. S = phi0/Kf + (alpha - phi0)(1 - alpha)/K = 0.0625
# Consolidation coefficient. c = k/(S + alpha^2 m) = 13.95348837
# Normal stress on top. q = 1
# Initial porepressure, resulting from instantaneous application of q, assuming corresponding instantaneous increase of porepressure (Note that this is calculated by MOOSE: we only need it for the analytical solution). p0 = alpha*m*q/(S + alpha^2 m) = 0.69767442
# Initial vertical displacement (down is positive), resulting from instantaneous application of q (Note this is calculated by MOOSE: we only need it for the analytical solution). uz0 = q*m*h*S/(S + alpha^2 m)
# Final vertical displacement (down in positive) (Note this is calculated by MOOSE: we only need it for the analytical solution). uzinf = q*m*h
#
# The solution for porepressure is
# P = 4*p0/\pi \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{2k-1} \cos ((2k-1)\pi z/(2h)) \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
# This series converges very slowly for ct/h^2 small, so in that domain
# P = p0 erf( (1-(z/h))/(2 \sqrt(ct/h^2)) )
#
# The degree of consolidation is defined as
# U = (uz - uz0)/(uzinf - uz0)
# where uz0 and uzinf are defined above, and
# uz = the vertical displacement of the top (down is positive)
# U = 1 - (8/\pi^2)\sum_{k=1}^{\infty} \frac{1}{(2k-1)^2} \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 10
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = 0
zmax = 10
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[../]
[./confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[../]
[./basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[../]
[./topdrained]
type = DirichletBC
variable = porepressure
value = 0
boundary = front
[../]
[./topload]
type = NeumannBC
variable = disp_z
value = -1
boundary = front
[../]
[]
[Modules]
[./FluidProperties]
[./the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 8.0
viscosity = 0.96
density0 = 1.0
[../]
[../]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
multiply_by_density = false
porepressure = porepressure
biot_coefficient = 0.6
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '2 3'
# bulk modulus is lambda + 2*mu/3 = 2 + 2*3/3 = 4
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./porosity]
type = PorousFlowPorosityConst # only the initial value of this is used
porosity = 0.1
[../]
[./biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
fluid_bulk_modulus = 8
solid_bulk_compliance = 0.25
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
use_displaced_mesh = false
[../]
[./p1]
type = PointValue
outputs = csv
point = '0 0 1'
variable = porepressure
use_displaced_mesh = false
[../]
[./p2]
type = PointValue
outputs = csv
point = '0 0 2'
variable = porepressure
use_displaced_mesh = false
[../]
[./p3]
type = PointValue
outputs = csv
point = '0 0 3'
variable = porepressure
use_displaced_mesh = false
[../]
[./p4]
type = PointValue
outputs = csv
point = '0 0 4'
variable = porepressure
use_displaced_mesh = false
[../]
[./p5]
type = PointValue
outputs = csv
point = '0 0 5'
variable = porepressure
use_displaced_mesh = false
[../]
[./p6]
type = PointValue
outputs = csv
point = '0 0 6'
variable = porepressure
use_displaced_mesh = false
[../]
[./p7]
type = PointValue
outputs = csv
point = '0 0 7'
variable = porepressure
use_displaced_mesh = false
[../]
[./p8]
type = PointValue
outputs = csv
point = '0 0 8'
variable = porepressure
use_displaced_mesh = false
[../]
[./p9]
type = PointValue
outputs = csv
point = '0 0 9'
variable = porepressure
use_displaced_mesh = false
[../]
[./p99]
type = PointValue
outputs = csv
point = '0 0 10'
variable = porepressure
use_displaced_mesh = false
[../]
[./zdisp]
type = PointValue
outputs = csv
point = '0 0 10'
variable = disp_z
use_displaced_mesh = false
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.5*t<0.1,0.5*t,0.1)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.0001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = terzaghi_basicthm
[./csv]
type = CSV
[../]
[]
modules/porous_flow/test/tests/poro_elasticity/terzaghi_constM.i
# Terzaghi's problem of consolodation of a drained medium
#
# A saturated soil sample sits in a bath of water.
# It is constrained on its sides, and bottom.
# Its sides and bottom are also impermeable.
# Initially it is unstressed.
# A normal stress, q, is applied to the soil's top.
# The soil then slowly compresses as water is squeezed
# out from the sample from its top (the top BC for
# the porepressure is porepressure = 0).
#
# See, for example. Section 2.2 of the online manuscript
# Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013
# but note that the "sigma" in that paper is the negative
# of the stress in TensorMechanics
#
# Here are the problem's parameters, and their values:
# Soil height. h = 10
# Soil's Lame lambda. la = 2
# Soil's Lame mu, which is also the Soil's shear modulus. mu = 3
# Soil bulk modulus. K = la + 2*mu/3 = 4
# Soil confined compressibility. m = 1/(K + 4mu/3) = 0.125
# Soil bulk compliance. 1/K = 0.25
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Soil initial storativity, which is the reciprocal of the initial Biot modulus. S = phi0/Kf + (alpha - phi0)(1 - alpha)/K = 0.0625
# Consolidation coefficient. c = k/(S + alpha^2 m) = 13.95348837
# Normal stress on top. q = 1
# Initial porepressure, resulting from instantaneous application of q, assuming corresponding instantaneous increase of porepressure (Note that this is calculated by MOOSE: we only need it for the analytical solution). p0 = alpha*m*q/(S + alpha^2 m) = 0.69767442
# Initial vertical displacement (down is positive), resulting from instantaneous application of q (Note this is calculated by MOOSE: we only need it for the analytical solution). uz0 = q*m*h*S/(S + alpha^2 m)
# Final vertical displacement (down in positive) (Note this is calculated by MOOSE: we only need it for the analytical solution). uzinf = q*m*h
#
# The solution for porepressure is
# P = 4*p0/\pi \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{2k-1} \cos ((2k-1)\pi z/(2h)) \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
# This series converges very slowly for ct/h^2 small, so in that domain
# P = p0 erf( (1-(z/h))/(2 \sqrt(ct/h^2)) )
#
# The degree of consolidation is defined as
# U = (uz - uz0)/(uzinf - uz0)
# where uz0 and uzinf are defined above, and
# uz = the vertical displacement of the top (down is positive)
# U = 1 - (8/\pi^2)\sum_{k=1}^{\infty} \frac{1}{(2k-1)^2} \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 10
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = 0
zmax = 10
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[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.8
alpha = 1
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[../]
[./confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[../]
[./basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[../]
[./topdrained]
type = DirichletBC
variable = porepressure
value = 0
boundary = front
[../]
[./topload]
type = NeumannBC
variable = disp_z
value = -1
boundary = front
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[../]
[./poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[../]
[./mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[../]
[./flux]
type = PorousFlowAdvectiveFlux
variable = porepressure
gravity = '0 0 0'
fluid_component = 0
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 0.96
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '2 3'
# bulk modulus is lambda + 2*mu/3 = 2 + 2*3/3 = 4
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./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]
type = PorousFlowPorosityHMBiotModulus
porosity_zero = 0.1
biot_coefficient = 0.6
solid_bulk = 4
constant_fluid_bulk_modulus = 8
constant_biot_modulus = 16
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0 # unimportant in this fully-saturated situation
phase = 0
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
use_displaced_mesh = false
[../]
[./p1]
type = PointValue
outputs = csv
point = '0 0 1'
variable = porepressure
use_displaced_mesh = false
[../]
[./p2]
type = PointValue
outputs = csv
point = '0 0 2'
variable = porepressure
use_displaced_mesh = false
[../]
[./p3]
type = PointValue
outputs = csv
point = '0 0 3'
variable = porepressure
use_displaced_mesh = false
[../]
[./p4]
type = PointValue
outputs = csv
point = '0 0 4'
variable = porepressure
use_displaced_mesh = false
[../]
[./p5]
type = PointValue
outputs = csv
point = '0 0 5'
variable = porepressure
use_displaced_mesh = false
[../]
[./p6]
type = PointValue
outputs = csv
point = '0 0 6'
variable = porepressure
use_displaced_mesh = false
[../]
[./p7]
type = PointValue
outputs = csv
point = '0 0 7'
variable = porepressure
use_displaced_mesh = false
[../]
[./p8]
type = PointValue
outputs = csv
point = '0 0 8'
variable = porepressure
use_displaced_mesh = false
[../]
[./p9]
type = PointValue
outputs = csv
point = '0 0 9'
variable = porepressure
use_displaced_mesh = false
[../]
[./p99]
type = PointValue
outputs = csv
point = '0 0 10'
variable = porepressure
use_displaced_mesh = false
[../]
[./zdisp]
type = PointValue
outputs = csv
point = '0 0 10'
variable = disp_z
use_displaced_mesh = false
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.5*t<0.1,0.5*t,0.1)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.0001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = terzaghi_constM
[./csv]
type = CSV
[../]
[]
modules/combined/test/tests/poro_mechanics/borehole_lowres.i
# Poroelastic response of a borehole.
#
# LOWRES VERSION: this version does not give perfect agreement with the analytical solution
#
# A fully-saturated medium contains a fluid with a homogeneous porepressure,
# but an anisitropic insitu stress. A infinitely-long borehole aligned with
# the $$z$$ axis is instanteously excavated. The borehole boundary is
# stress-free and allowed to freely drain. This problem is analysed using
# plane-strain conditions (no $$z$$ displacement).
#
# The solution in Laplace space is found in E Detournay and AHD Cheng "Poroelastic response of a borehole in a non-hydrostatic stress field". International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 25 (1988) 171-182. In the small-time limit, the Laplace transforms may be performed. There is one typo in the paper. Equation (A4)'s final term should be -(a/r)\sqrt(4ct/(a^2\pi)), and not +(a/r)\sqrt(4ct/(a^2\pi)).
#
# Because realistic parameters are chosen (below),
# the residual for porepressure is much smaller than
# the residuals for the displacements. Therefore the
# scaling parameter is chosen. Also note that the
# insitu stresses are effective stresses, not total
# stresses, but the solution in the above paper is
# expressed in terms of total stresses.
#
# Here are the problem's parameters, and their values:
# Borehole radius. a = 1
# Rock's Lame lambda. la = 0.5E9
# Rock's Lame mu, which is also the Rock's shear modulus. mu = G = 1.5E9
# Rock bulk modulus. K = la + 2*mu/3 = 1.5E9
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.125
# Rock bulk compliance. 1/K = 0.66666666E-9
# Fluid bulk modulus. Kf = 0.7171315E9
# Fluid bulk compliance. 1/Kf = 1.39444444E-9
# Rock initial porosity. phi0 = 0.3
# Biot coefficient. alpha = 0.65
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 2E9
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.345E9
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.2364
# Skempton coefficient. B = alpha*M/Ku = 0.554
# Fluid mobility (rock permeability/fluid viscosity). k = 1E-12
[Mesh]
type = FileMesh
file = borehole_lowres_input.e
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
block = 1
[]
[GlobalParams]
volumetric_locking_correction=true
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
scaling = 1E9 # Notice the scaling, to make porepressure's kernels roughly of same magnitude as disp's kernels
[../]
[]
[ICs]
[./initial_p]
type = ConstantIC
variable = porepressure
value = 1E6
[../]
[]
[BCs]
[./fixed_outer_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = outer
[../]
[./fixed_outer_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = outer
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'zmin zmax'
[../]
[./borehole_wall]
type = DirichletBC
variable = porepressure
value = 0
boundary = bh_wall
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./tot_yy]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./tot_yy]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_yy
function = 'stress_yy-0.65*porepressure'
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PoroMechanicsCoupling
variable = disp_x
component = 0
[../]
[./poro_y]
type = PoroMechanicsCoupling
variable = disp_y
component = 1
[../]
[./poro_z]
type = PoroMechanicsCoupling
variable = disp_z
component = 2
[../]
[./poro_timederiv]
type = PoroFullSatTimeDerivative
variable = porepressure
[../]
[./darcy_flow]
type = CoefDiffusion
variable = porepressure
coef = 1E-12
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5E9 1.5E9'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*1.5/3 = 1.5E9
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeFiniteStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '-1.35E6 0 0 0 -3.35E6 0 0 0 0' # remember this is the effective stress
eigenstrain_name = ini_stress
[../]
[./no_plasticity]
type = ComputeFiniteStrainElasticStress
[../]
[./poro_material]
type = PoroFullSatMaterial
porosity0 = 0.3
biot_coefficient = 0.65
solid_bulk_compliance = 0.6666666666667E-9
fluid_bulk_compliance = 1.3944444444444E-9
constant_porosity = false
[../]
[]
[Postprocessors]
[./p00]
type = PointValue
variable = porepressure
point = '1.00 0 0'
outputs = csv_p
[../]
[./p01]
type = PointValue
variable = porepressure
point = '1.01 0 0'
outputs = csv_p
[../]
[./p02]
type = PointValue
variable = porepressure
point = '1.02 0 0'
outputs = csv_p
[../]
[./p03]
type = PointValue
variable = porepressure
point = '1.03 0 0'
outputs = csv_p
[../]
[./p04]
type = PointValue
variable = porepressure
point = '1.04 0 0'
outputs = csv_p
[../]
[./p05]
type = PointValue
variable = porepressure
point = '1.05 0 0'
outputs = csv_p
[../]
[./p06]
type = PointValue
variable = porepressure
point = '1.06 0 0'
outputs = csv_p
[../]
[./p07]
type = PointValue
variable = porepressure
point = '1.07 0 0'
outputs = csv_p
[../]
[./p08]
type = PointValue
variable = porepressure
point = '1.08 0 0'
outputs = csv_p
[../]
[./p09]
type = PointValue
variable = porepressure
point = '1.09 0 0'
outputs = csv_p
[../]
[./p10]
type = PointValue
variable = porepressure
point = '1.10 0 0'
outputs = csv_p
[../]
[./p11]
type = PointValue
variable = porepressure
point = '1.11 0 0'
outputs = csv_p
[../]
[./p12]
type = PointValue
variable = porepressure
point = '1.12 0 0'
outputs = csv_p
[../]
[./p13]
type = PointValue
variable = porepressure
point = '1.13 0 0'
outputs = csv_p
[../]
[./p14]
type = PointValue
variable = porepressure
point = '1.14 0 0'
outputs = csv_p
[../]
[./p15]
type = PointValue
variable = porepressure
point = '1.15 0 0'
outputs = csv_p
[../]
[./p16]
type = PointValue
variable = porepressure
point = '1.16 0 0'
outputs = csv_p
[../]
[./p17]
type = PointValue
variable = porepressure
point = '1.17 0 0'
outputs = csv_p
[../]
[./p18]
type = PointValue
variable = porepressure
point = '1.18 0 0'
outputs = csv_p
[../]
[./p19]
type = PointValue
variable = porepressure
point = '1.19 0 0'
outputs = csv_p
[../]
[./p20]
type = PointValue
variable = porepressure
point = '1.20 0 0'
outputs = csv_p
[../]
[./p21]
type = PointValue
variable = porepressure
point = '1.21 0 0'
outputs = csv_p
[../]
[./p22]
type = PointValue
variable = porepressure
point = '1.22 0 0'
outputs = csv_p
[../]
[./p23]
type = PointValue
variable = porepressure
point = '1.23 0 0'
outputs = csv_p
[../]
[./p24]
type = PointValue
variable = porepressure
point = '1.24 0 0'
outputs = csv_p
[../]
[./p25]
type = PointValue
variable = porepressure
point = '1.25 0 0'
outputs = csv_p
[../]
[./s00]
type = PointValue
variable = disp_x
point = '1.00 0 0'
outputs = csv_s
[../]
[./s01]
type = PointValue
variable = disp_x
point = '1.01 0 0'
outputs = csv_s
[../]
[./s02]
type = PointValue
variable = disp_x
point = '1.02 0 0'
outputs = csv_s
[../]
[./s03]
type = PointValue
variable = disp_x
point = '1.03 0 0'
outputs = csv_s
[../]
[./s04]
type = PointValue
variable = disp_x
point = '1.04 0 0'
outputs = csv_s
[../]
[./s05]
type = PointValue
variable = disp_x
point = '1.05 0 0'
outputs = csv_s
[../]
[./s06]
type = PointValue
variable = disp_x
point = '1.06 0 0'
outputs = csv_s
[../]
[./s07]
type = PointValue
variable = disp_x
point = '1.07 0 0'
outputs = csv_s
[../]
[./s08]
type = PointValue
variable = disp_x
point = '1.08 0 0'
outputs = csv_s
[../]
[./s09]
type = PointValue
variable = disp_x
point = '1.09 0 0'
outputs = csv_s
[../]
[./s10]
type = PointValue
variable = disp_x
point = '1.10 0 0'
outputs = csv_s
[../]
[./s11]
type = PointValue
variable = disp_x
point = '1.11 0 0'
outputs = csv_s
[../]
[./s12]
type = PointValue
variable = disp_x
point = '1.12 0 0'
outputs = csv_s
[../]
[./s13]
type = PointValue
variable = disp_x
point = '1.13 0 0'
outputs = csv_s
[../]
[./s14]
type = PointValue
variable = disp_x
point = '1.14 0 0'
outputs = csv_s
[../]
[./s15]
type = PointValue
variable = disp_x
point = '1.15 0 0'
outputs = csv_s
[../]
[./s16]
type = PointValue
variable = disp_x
point = '1.16 0 0'
outputs = csv_s
[../]
[./s17]
type = PointValue
variable = disp_x
point = '1.17 0 0'
outputs = csv_s
[../]
[./s18]
type = PointValue
variable = disp_x
point = '1.18 0 0'
outputs = csv_s
[../]
[./s19]
type = PointValue
variable = disp_x
point = '1.19 0 0'
outputs = csv_s
[../]
[./s20]
type = PointValue
variable = disp_x
point = '1.20 0 0'
outputs = csv_s
[../]
[./s21]
type = PointValue
variable = disp_x
point = '1.21 0 0'
outputs = csv_s
[../]
[./s22]
type = PointValue
variable = disp_x
point = '1.22 0 0'
outputs = csv_s
[../]
[./s23]
type = PointValue
variable = disp_x
point = '1.23 0 0'
outputs = csv_s
[../]
[./s24]
type = PointValue
variable = disp_x
point = '1.24 0 0'
outputs = csv_s
[../]
[./s25]
type = PointValue
variable = disp_x
point = '1.25 0 0'
outputs = csv_s
[../]
[./t00]
type = PointValue
variable = tot_yy
point = '1.00 0 0'
outputs = csv_t
[../]
[./t01]
type = PointValue
variable = tot_yy
point = '1.01 0 0'
outputs = csv_t
[../]
[./t02]
type = PointValue
variable = tot_yy
point = '1.02 0 0'
outputs = csv_t
[../]
[./t03]
type = PointValue
variable = tot_yy
point = '1.03 0 0'
outputs = csv_t
[../]
[./t04]
type = PointValue
variable = tot_yy
point = '1.04 0 0'
outputs = csv_t
[../]
[./t05]
type = PointValue
variable = tot_yy
point = '1.05 0 0'
outputs = csv_t
[../]
[./t06]
type = PointValue
variable = tot_yy
point = '1.06 0 0'
outputs = csv_t
[../]
[./t07]
type = PointValue
variable = tot_yy
point = '1.07 0 0'
outputs = csv_t
[../]
[./t08]
type = PointValue
variable = tot_yy
point = '1.08 0 0'
outputs = csv_t
[../]
[./t09]
type = PointValue
variable = tot_yy
point = '1.09 0 0'
outputs = csv_t
[../]
[./t10]
type = PointValue
variable = tot_yy
point = '1.10 0 0'
outputs = csv_t
[../]
[./t11]
type = PointValue
variable = tot_yy
point = '1.11 0 0'
outputs = csv_t
[../]
[./t12]
type = PointValue
variable = tot_yy
point = '1.12 0 0'
outputs = csv_t
[../]
[./t13]
type = PointValue
variable = tot_yy
point = '1.13 0 0'
outputs = csv_t
[../]
[./t14]
type = PointValue
variable = tot_yy
point = '1.14 0 0'
outputs = csv_t
[../]
[./t15]
type = PointValue
variable = tot_yy
point = '1.15 0 0'
outputs = csv_t
[../]
[./t16]
type = PointValue
variable = tot_yy
point = '1.16 0 0'
outputs = csv_t
[../]
[./t17]
type = PointValue
variable = tot_yy
point = '1.17 0 0'
outputs = csv_t
[../]
[./t18]
type = PointValue
variable = tot_yy
point = '1.18 0 0'
outputs = csv_t
[../]
[./t19]
type = PointValue
variable = tot_yy
point = '1.19 0 0'
outputs = csv_t
[../]
[./t20]
type = PointValue
variable = tot_yy
point = '1.20 0 0'
outputs = csv_t
[../]
[./t21]
type = PointValue
variable = tot_yy
point = '1.21 0 0'
outputs = csv_t
[../]
[./t22]
type = PointValue
variable = tot_yy
point = '1.22 0 0'
outputs = csv_t
[../]
[./t23]
type = PointValue
variable = tot_yy
point = '1.23 0 0'
outputs = csv_t
[../]
[./t24]
type = PointValue
variable = tot_yy
point = '1.24 0 0'
outputs = csv_t
[../]
[./t25]
type = PointValue
variable = tot_yy
point = '1.25 0 0'
outputs = csv_t
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = 2*t
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options = '-snes_monitor -snes_linesearch_monitor'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'gmres asm 1E0 1E-10 200 500 lu NONZERO'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.3
dt = 0.3
#[./TimeStepper]
# type = PostprocessorDT
# postprocessor = dt
# dt = 0.003
#[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = borehole_lowres
exodus = true
sync_times = '0.003 0.3'
[./csv_p]
file_base = borehole_lowres_p
type = CSV
[../]
[./csv_s]
file_base = borehole_lowres_s
type = CSV
[../]
[./csv_t]
file_base = borehole_lowres_t
type = CSV
[../]
[]
modules/porous_flow/test/tests/poro_elasticity/mandel_constM.i
# Mandel's problem of consolodation of a drained medium
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[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.8
alpha = 1e-5
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[../]
[./roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[../]
[./xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[../]
[./top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[../]
[]
[Functions]
[./top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./tot_force]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./tot_force]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
function = '-stress_yy+0.6*porepressure'
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[../]
[./poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[../]
[./mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[../]
[./flux]
type = PorousFlowAdvectiveFlux
variable = porepressure
gravity = '0 0 0'
fluid_component = 0
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 1
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./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]
type = PorousFlowPorosityHMBiotModulus
porosity_zero = 0.1
biot_coefficient = 0.6
solid_bulk = 1
constant_fluid_bulk_modulus = 8
constant_biot_modulus = 4.7058823529
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0 # unimportant in this fully-saturated situation
phase = 0
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[../]
[./p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[../]
[./p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[../]
[./p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[../]
[./p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[../]
[./p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[../]
[./p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[../]
[./p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[../]
[./p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[../]
[./p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[../]
[./p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[../]
[./xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[../]
[./ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[../]
[./total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel_constM
[./csv]
interval = 3
type = CSV
[../]
[]
modules/level_set/examples/rotating_circle/circle_rotate_supg.i
[Mesh]
type = GeneratedMesh
dim = 2
xmin = -1
xmax = 1
ymin = -1
ymax = 1
nx = 32
ny = 32
uniform_refine = 2
[]
[AuxVariables]
[./vel_x]
[../]
[./vel_y]
[../]
[]
[AuxKernels]
[./vel_x]
type = FunctionAux
function = 4*y
variable = vel_x
execute_on = initial
[../]
[./vel_y]
type = FunctionAux
function = -4*x
variable = vel_y
execute_on = initial
[../]
[]
[Variables]
[./phi]
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.03
center = '0 0.5 0'
radius = 0.15
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[./advection_supg]
type = LevelSetAdvectionSUPG
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[./time_supg]
type = LevelSetTimeDerivativeSUPG
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[../]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 1.570796
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[../]
[]
[Outputs]
csv = true
exodus = true
[]
modules/level_set/examples/vortex/vortex_reinit.i
[Mesh]
type = GeneratedMesh
dim = 2
xmax = 1
ymax = 1
nx = 16
ny = 16
uniform_refine = 2
elem_type = QUAD9
second_order = true
[]
[AuxVariables]
[vel_x]
family = LAGRANGE
[]
[vel_y]
family = LAGRANGE
[]
[]
[AuxKernels]
[vel_x]
type = FunctionAux
function = vel_x
variable = vel_x
execute_on = 'initial timestep_begin'
[]
[vel_y]
type = FunctionAux
function = vel_y
variable = vel_y
execute_on = 'initial timestep_begin'
[]
[]
[Variables]
[phi]
family = LAGRANGE
[]
[]
[Functions]
[phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.03
center = '0.5 0.75 0'
radius = 0.15
[]
[vel_x]
type = LevelSetOlssonVortex
component = x
reverse_time = 2
[]
[vel_y]
type = LevelSetOlssonVortex
component = y
reverse_time = 2
[]
[]
[ICs]
[phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[]
[]
[Kernels]
[time]
type = TimeDerivative
variable = phi
[]
[advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[]
[advection_supg]
type = LevelSetAdvectionSUPG
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[]
[time_supg]
type = LevelSetTimeDerivativeSUPG
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[]
[]
[Postprocessors]
[area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[]
[cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial timestep_end'
[]
[]
[Problem]
type = LevelSetProblem
[]
[Preconditioning/smp]
type = SMP
full = true
[]
[Executioner]
type = Transient
solve_type = NEWTON
start_time = 0
end_time = 2
scheme = crank-nicolson
[TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[]
[]
[MultiApps]
[reinit]
type = LevelSetReinitializationMultiApp
input_files = 'vortex_reinit_sub.i'
execute_on = TIMESTEP_END
[]
[]
[Transfers]
[to_sub]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[]
[to_sub_init]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi_0
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[]
[from_sub]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi
direction = from_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[]
[]
[Outputs]
csv = true
exodus = true
[]
modules/combined/test/tests/poro_mechanics/terzaghi.i
# Terzaghi's problem of consolodation of a drained medium
#
# A saturated soil sample sits in a bath of water.
# It is constrained on its sides, and bottom.
# Its sides and bottom are also impermeable.
# Initially it is unstressed.
# A normal stress, q, is applied to the soil's top.
# The soil then slowly compresses as water is squeezed
# out from the sample from its top (the top BC for
# the porepressure is porepressure = 0).
#
# See, for example. Section 2.2 of the online manuscript
# Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013
# but note that the "sigma" in that paper is the negative
# of the stress in TensorMechanics
#
# Here are the problem's parameters, and their values:
# Soil height. h = 10
# Soil's Lame lambda. la = 2
# Soil's Lame mu, which is also the Soil's shear modulus. mu = 3
# Soil bulk modulus. K = la + 2*mu/3 = 4
# Soil confined compressibility. m = 1/(K + 4mu/3) = 0.125
# Soil bulk compliance. 1/K = 0.25
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Soil initial storativity, which is the reciprocal of the initial Biot modulus. S = phi0/Kf + (alpha - phi0)(1 - alpha)/K = 0.0625
# Consolidation coefficient. c = k/(S + alpha^2 m) = 13.95348837
# Normal stress on top. q = 1
# Initial porepressure, resulting from instantaneous application of q, assuming corresponding instantaneous increase of porepressure (Note that this is calculated by MOOSE: we only need it for the analytical solution). p0 = alpha*m*q/(S + alpha^2 m) = 0.69767442
# Initial vertical displacement (down is positive), resulting from instantaneous application of q (Note this is calculated by MOOSE: we only need it for the analytical solution). uz0 = q*m*h*S/(S + alpha^2 m)
# Final vertical displacement (down in positive) (Note this is calculated by MOOSE: we only need it for the analytical solution). uzinf = q*m*h
#
# The solution for porepressure is
# P = 4*p0/\pi \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{2k-1} \cos ((2k-1)\pi z/(2h)) \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
# This series converges very slowly for ct/h^2 small, so in that domain
# P = p0 erf( (1-(z/h))/(2 \sqrt(ct/h^2)) )
#
# The degree of consolidation is defined as
# U = (uz - uz0)/(uzinf - uz0)
# where uz0 and uzinf are defined above, and
# uz = the vertical displacement of the top (down is positive)
# U = 1 - (8/\pi^2)\sum_{k=1}^{\infty} \frac{1}{(2k-1)^2} \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 10
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = 0
zmax = 10
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[../]
[./confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[../]
[./basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[../]
[./topdrained]
type = DirichletBC
variable = porepressure
value = 0
boundary = front
[../]
[./topload]
type = NeumannBC
variable = disp_z
value = -1
boundary = front
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PoroMechanicsCoupling
variable = disp_x
component = 0
[../]
[./poro_y]
type = PoroMechanicsCoupling
variable = disp_y
component = 1
[../]
[./poro_z]
type = PoroMechanicsCoupling
variable = disp_z
component = 2
[../]
[./poro_timederiv]
type = PoroFullSatTimeDerivative
variable = porepressure
[../]
[./darcy_flow]
type = CoefDiffusion
variable = porepressure
coef = 1.5
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '2 3'
# bulk modulus is lambda + 2*mu/3 = 2 + 2*3/3 = 4
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./poro_material]
type = PoroFullSatMaterial
porosity0 = 0.1
biot_coefficient = 0.6
solid_bulk_compliance = 0.25
fluid_bulk_compliance = 0.125
constant_porosity = true
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
[../]
[./p1]
type = PointValue
outputs = csv
point = '0 0 1'
variable = porepressure
[../]
[./p2]
type = PointValue
outputs = csv
point = '0 0 2'
variable = porepressure
[../]
[./p3]
type = PointValue
outputs = csv
point = '0 0 3'
variable = porepressure
[../]
[./p4]
type = PointValue
outputs = csv
point = '0 0 4'
variable = porepressure
[../]
[./p5]
type = PointValue
outputs = csv
point = '0 0 5'
variable = porepressure
[../]
[./p6]
type = PointValue
outputs = csv
point = '0 0 6'
variable = porepressure
[../]
[./p7]
type = PointValue
outputs = csv
point = '0 0 7'
variable = porepressure
[../]
[./p8]
type = PointValue
outputs = csv
point = '0 0 8'
variable = porepressure
[../]
[./p9]
type = PointValue
outputs = csv
point = '0 0 9'
variable = porepressure
[../]
[./p99]
type = PointValue
outputs = csv
point = '0 0 10'
variable = porepressure
[../]
[./zdisp]
type = PointValue
outputs = csv
point = '0 0 10'
variable = disp_z
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.5*t<0.1,0.5*t,0.1)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-14 1E-10 10000'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.0001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = terzaghi
[./csv]
type = CSV
[../]
[]
modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated_volume.i
# Mandel's problem of consolodation of a drained medium
# Using the FullySaturatedDarcyBase and FullySaturatedFullySaturatedMassTimeDerivative kernels
# with multiply_by_density = false, so that this problem becomes linear
# Note the use of consistent_with_displaced_mesh = false in the calculation of volumetric strain
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[../]
[./roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[../]
[./xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[../]
[./top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[../]
[]
[Functions]
[./top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./tot_force]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./tot_force]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
function = '-stress_yy+0.6*porepressure'
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[../]
[./mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
biot_coefficient = 0.6
multiply_by_density = false
coupling_type = HydroMechanical
variable = porepressure
[../]
[./flux]
type = PorousFlowFullySaturatedDarcyBase
multiply_by_density = false
variable = porepressure
gravity = '0 0 0'
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 1
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./eff_fluid_pressure_qp]
type = PorousFlowEffectiveFluidPressure
[../]
[./vol_strain]
type = PorousFlowVolumetricStrain
consistent_with_displaced_mesh = false
[../]
[./ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = porepressure
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid_qp]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[../]
[./porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[../]
[./biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
solid_bulk_compliance = 1
fluid_bulk_modulus = 8
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[../]
[./p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[../]
[./p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[../]
[./p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[../]
[./p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[../]
[./p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[../]
[./p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[../]
[./p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[../]
[./p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[../]
[./p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[../]
[./p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[../]
[./xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[../]
[./ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[../]
[./total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel_fully_saturated_volume
[./csv]
interval = 3
type = CSV
[../]
[]
modules/level_set/examples/circle/circle_16.i
[Mesh]
type = GeneratedMesh
dim = 2
xmin = 0
xmax = 1
ymin = 0
ymax = 1
nx = 16
ny = 16
[]
[Variables]
[./phi]
[../]
[]
[AuxVariables]
[./vel_x]
initial_condition = 3
[../]
[./vel_y]
initial_condition = 3
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.05
center = '0.5 0.5 0'
radius = 0.15
[../]
[]
[BCs]
[./Periodic]
[./all]
variable = phi
auto_direction = 'x y'
[../]
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 1
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[../]
[]
[Outputs]
csv = true
exodus = true
[]
modules/level_set/examples/rotating_circle/circle_rotate_master.i
[Mesh]
type = GeneratedMesh
dim = 2
xmin = -1
xmax = 1
ymin = -1
ymax = 1
nx = 32
ny = 32
uniform_refine = 2
[]
[AuxVariables]
[./vel_x]
[../]
[./vel_y]
[../]
[]
[AuxKernels]
[./vel_x]
type = FunctionAux
function = 4*y
variable = vel_x
execute_on = initial
[../]
[./vel_y]
type = FunctionAux
function = -4*x
variable = vel_y
execute_on = initial
[../]
[]
[Variables]
[./phi]
[../]
[]
[BCs]
[./all]
type = DirichletBC
variable = phi
boundary = 'top bottom left right'
value = 0
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.03
center = '0 0.5 0'
radius = 0.15
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[../]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial timestep_end'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 1.570796
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[../]
[]
[MultiApps]
[./reinit]
type = LevelSetReinitializationMultiApp
input_files = 'circle_rotate_sub.i'
execute_on = 'timestep_end'
[../]
[]
[Transfers]
[./to_sub]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[./to_sub_init]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi_0
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[./from_sub]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi
direction = from_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[]
[Outputs]
csv = true
exodus = true
[]
modules/porous_flow/test/tests/poro_elasticity/mandel_basicthm.i
# using a BasicTHM Action
#
# Mandel's problem of consolodation of a drained medium
# Using the FullySaturatedDarcyBase and FullySaturatedFullySaturatedMassTimeDerivative kernels
# with multiply_by_density = false, so that this problem becomes linear
# Note the use of consistent_with_displaced_mesh = false in the calculation of volumetric strain
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[../]
[./roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[../]
[./xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[../]
[./top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[../]
[]
[Functions]
[./top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[../]
[]
[AuxVariables]
[./tot_force]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./tot_force]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
function = '-stress_yy+0.6*porepressure'
[../]
[]
[Modules]
[./FluidProperties]
[./the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 8.0
viscosity = 1.0
density0 = 1.0
[../]
[../]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
multiply_by_density = false
porepressure = porepressure
biot_coefficient = 0.6
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[../]
[./biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
solid_bulk_compliance = 1
fluid_bulk_modulus = 8
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[../]
[./p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[../]
[./p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[../]
[./p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[../]
[./p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[../]
[./p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[../]
[./p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[../]
[./p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[../]
[./p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[../]
[./p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[../]
[./p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[../]
[./xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[../]
[./ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[../]
[./total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel_basicthm
[./csv]
interval = 3
type = CSV
[../]
[]
test/tests/time_steppers/postprocessor_dt/postprocessor_dt.i
[Mesh]
type = GeneratedMesh
dim = 2
xmin = -1
xmax = 1
ymin = -1
ymax = 1
nx = 2
ny = 2
elem_type = QUAD9
[]
[Functions]
[./exact_fn]
type = ParsedFunction
value = t*t*(x*x+y*y)
[../]
[./forcing_fn]
type = ParsedFunction
value = 2*t*(x*x+y*y)-4*t*t
[../]
[]
[Variables]
[./u]
family = LAGRANGE
order = SECOND
[../]
[]
[ICs]
[./u_var]
type = FunctionIC
variable = u
function = exact_fn
[../]
[]
[Kernels]
[./td]
type = TimeDerivative
variable = u
[../]
[./diff]
type = Diffusion
variable = u
[../]
[./ffn]
type = BodyForce
variable = u
function = forcing_fn
[../]
[]
[BCs]
[./all]
type = FunctionDirichletBC
variable = u
boundary = 'left right top bottom'
function = exact_fn
[../]
[]
[Postprocessors]
[./l2_error]
type = ElementL2Error
variable = u
function = exact_fn
execute_on = 'initial timestep_end'
[../]
# Just use some postprocessor that gives values good enough for time stepping ;-)
[./dt]
type = ElementAverageValue
variable = u
execute_on = 'initial timestep_end'
[../]
[]
[Executioner]
type = Transient
scheme = 'crank-nicolson'
start_time = 1.0
num_steps = 2
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
[../]
[]
[Outputs]
exodus = true
[]
modules/porous_flow/test/tests/poro_elasticity/mandel.i
# Mandel's problem of consolodation of a drained medium
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
#
# FINAL NOTE: The above solution assumes constant Biot Modulus.
# In porous_flow this is not true. Therefore the solution is
# a little different than in the paper. This test was therefore
# validated against MOOSE's poromechanics, which can choose either
# a constant Biot Modulus (which has been shown to agree with
# the analytic solution), or a non-constant Biot Modulus (which
# gives the same results as porous_flow).
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[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.8
alpha = 1e-5
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[../]
[./roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[../]
[./xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[../]
[./top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[../]
[]
[Functions]
[./top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./tot_force]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./tot_force]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
function = '-stress_yy+0.6*porepressure'
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[../]
[./poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[../]
[./mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[../]
[./flux]
type = PorousFlowAdvectiveFlux
variable = porepressure
gravity = '0 0 0'
fluid_component = 0
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 1
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./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]
type = PorousFlowPorosity
fluid = true
mechanical = true
ensure_positive = false
porosity_zero = 0.1
biot_coefficient = 0.6
solid_bulk = 1
[../]
[./permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0 # unimportant in this fully-saturated situation
phase = 0
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[../]
[./p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[../]
[./p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[../]
[./p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[../]
[./p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[../]
[./p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[../]
[./p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[../]
[./p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[../]
[./p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[../]
[./p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[../]
[./p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[../]
[./xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[../]
[./ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[../]
[./total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel
[./csv]
interval = 3
type = CSV
[../]
[]
modules/combined/test/tests/poro_mechanics/mandel.i
# Mandel's problem of consolodation of a drained medium
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./porepressure]
[../]
[]
[BCs]
[./roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[../]
[./roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[../]
[./plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[../]
[./xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[../]
[./top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[../]
[]
[Functions]
[./top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./tot_force]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./tot_force]
type = ParsedAux
args = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
function = '-stress_yy+0.6*porepressure'
[../]
[]
[Kernels]
[./grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[../]
[./grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[../]
[./grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[../]
[./poro_x]
type = PoroMechanicsCoupling
variable = disp_x
component = 0
[../]
[./poro_y]
type = PoroMechanicsCoupling
variable = disp_y
component = 1
[../]
[./poro_z]
type = PoroMechanicsCoupling
variable = disp_z
component = 2
[../]
[./poro_timederiv]
type = PoroFullSatTimeDerivative
variable = porepressure
[../]
[./darcy_flow]
type = CoefDiffusion
variable = porepressure
coef = 1.5
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[../]
[./strain]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./poro_material]
type = PoroFullSatMaterial
porosity0 = 0.1
biot_coefficient = 0.6
solid_bulk_compliance = 1
fluid_bulk_compliance = 0.125
constant_porosity = true
[../]
[]
[Postprocessors]
[./p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[../]
[./p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[../]
[./p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[../]
[./p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[../]
[./p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[../]
[./p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[../]
[./p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[../]
[./p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[../]
[./p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[../]
[./p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[../]
[./p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[../]
[./xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[../]
[./ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[../]
[./total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[../]
[./dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-14 1E-10 10000'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[./TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel
[./csv]
interval = 3
type = CSV
[../]
[]
modules/level_set/examples/vortex/vortex.i
[Mesh]
type = GeneratedMesh
dim = 2
xmax = 1
ymax = 1
nx = 16
ny = 16
uniform_refine = 2
elem_type = QUAD9
[]
[AuxVariables]
[./vel_x]
family = LAGRANGE
order = FIRST
[../]
[./vel_y]
family = LAGRANGE
order = FIRST
[../]
[]
[AuxKernels]
[./vel_x]
type = FunctionAux
function = vel_x
variable = vel_x
execute_on = 'initial timestep_begin'
[../]
[./vel_y]
type = FunctionAux
function = vel_y
variable = vel_y
execute_on = 'initial timestep_begin'
[../]
[]
[Variables]
[./phi]
family = LAGRANGE
order = FIRST
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.01184
center = '0.5 0.75 0'
radius = 0.15
[../]
[./vel_x]
type = LevelSetOlssonVortex
component = x
reverse_time = 2
[../]
[./vel_y]
type = LevelSetOlssonVortex
component = y
reverse_time = 2
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[../]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial timestep_end'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 2
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[../]
[]
[Outputs]
csv = true
exodus = true
[]