HM_5km_SF.i

## Steady State detection - No injection/production
## Same elasticity tensor for the 3 layers
## Constant porosity and permeability
## Changed the bulk modulus of water from 2 to 6e9)

density_overburden = 2420

[Mesh]
  allow_renumbering = true
  [main]
    type = FileMeshGenerator
    file = 5km_mesh.msh
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  fp = the_simple_fluid
  displacements = 'disp_x disp_y disp_z'
[]

[Variables]
  [porepressure]
  []
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables]
  [stress_zz]
    family = MONOMIAL
    order = SECOND
  []
[]

[AuxKernels]
  [stress_zz]
    type = RankTwoAux
    variable = stress_zz
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  []
[]

[ICs]
  [Pressure]
    type = FunctionIC
    variable = porepressure
    function = 'hydrostaticP'
  []
  [stress]
    type = FunctionIC
    variable = stress_zz
    function = 'lithost_zz_BC'
  []
[]

[BCs]
  [Pore_Pressure_BC]
    type = FunctionDirichletBC
    variable = porepressure
    ##boundary = 'Top Bottom'
    ## boundary = 'Top Bottom Front Left'
    ## boundary = 'Top Bottom Front Back'
    ## boundary = 'Top Bottom Left Back Right Front'
    boundary = 'Bottom Left Back Right Front'
    function = 'hydrostaticP'
  []
  [Disp_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'Back Front'
    value = 0
  []
  [Disp_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'Left Right'
    value = 0
  []
  [Disp_z]
    type = DirichletBC
    variable = disp_z
    boundary = 'Bottom'
    value = 0
  []
  [litho_stress_zz]
    type = FunctionNeumannBC
    boundary = 'Top'
    variable = disp_z
    function = 'lithost_zz_BC' #Lithostatic pressure
  []
[]

[Functions]
  [hydrostaticP]
    type = ParsedFunction
    expression = '1000*9.8*(-z)+1e5'
  []
  [dts]
    type = PiecewiseLinear
    x = '0 3600 86400 864000 15768000
          15768001 15854400 16632000 31536000'
    y = '1 1800 43200 172800 864000
         3600 43200 86400 864000'
  []
  [lithost_zz_BC]
    type = ParsedFunction
    expression = '-${density_overburden}*9.81*(-z)'
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 6e9 ## Should be 2e9
    viscosity = 0.00026383
    thermal_expansion = 0
  []
[]

[Kernels]
  ##### Hydraulics #####
  [mass_dot]
    type = PorousFlowFullySaturatedMassTimeDerivative
    variable = porepressure
  []
  [advection]
    type = PorousFlowFullySaturatedAdvectiveFlux
    fluid_component = 0
    variable = porepressure
    gravity = '0 0 -9.81'
  []
  ####### Mechanics ########
  [gravity]
    type = Gravity
    use_displaced_mesh = false
    variable = disp_z
    value = -9.81 # MPa
    density = ${density_overburden}
  []
  [vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    variable = porepressure
  []
  [grad_stress_x]
    type = StressDivergenceTensors
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  []
  [grad_stress_y]
    type = StressDivergenceTensors
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  []
  [grad_stress_z]
    type = StressDivergenceTensors
    variable = disp_z
    use_displaced_mesh = false
    component = 2
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_z
    use_displaced_mesh = false
    component = 2
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_x disp_y disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [fluid_mass_prod]
    type = PorousFlowSumQuantity
  []
[]

[Materials]
  [undrained_density]
    type = PorousFlowTotalGravitationalDensityFullySaturatedFromPorosity
    rho_s = 3000 ## (por=15%, pf=1000, pb=2700)
    output_properties = density
    outputs = 'exodus'
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 1
    fluid_bulk_modulus = 2e9 #Keep it consistant with the FluidProperties
  []
  [saturation_calculator]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [water]
    type = PorousFlowSingleComponentFluid
    phase = 0
  []
  [porosity_const]
    type = PorousFlowPorosityConst
    porosity = 0.15
  []
  [permeability_res]
    type = PorousFlowPermeabilityConst
    permeability = '2.5e-14 0 0  0 2.5e-14 0  0 0 2.5e-14'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [eff_fluid_pressure]
    # Used by variant porosity
    type = PorousFlowEffectiveFluidPressure
  []
  [temperature]
    type = PorousFlowTemperature
    # temperature = temperature
  []
  ######## Mechanics ########
  [elasticity_tensor_res]
    type = ComputeIsotropicElasticityTensor
    bulk_modulus = 55e9
    shear_modulus = 32.59e9
  []
  [strain]
    type = ComputeSmallStrain
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[Preconditioning]
  active = 'HM'
  [smp]
    type = SMP
    full = true
  []
  [superlu]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'gmres lu superlu_dist'
  []
  [HM]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_hypre_type -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_agg_nl  -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_max_levels  -pc_hypre_boomeramg_coarsen_type   -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_P_max -pc_hypre_boomeramg_truncfactor'
    petsc_options_value = 'hypre    boomeramg 0.7 4 5 25 HMIS  ext+i 2 0.3'
  []
[]

[Executioner]
  type = Transient
  end_time = 31536000 #1 year (365 days)
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
  dtmin = 1
  nl_rel_tol = 1e-8
  nl_rel_step_tol = 1e-8
  l_tol = 1e-8
  automatic_scaling = true
  compute_scaling_once = true
  solve_type = NEWTON
  start_time = 0
  steady_state_detection = true
  steady_state_tolerance = 1e-09
[]

[Outputs]
  console = true
  csv = false
  exodus = true
  print_linear_residuals = false
[]

[Debug]
  show_var_residual_norms = true
[]