The calculation results differ between using and not using the PorousFlowUnsaturated action

[PEI0214]

Hi everyone,
I am trying to use the PorousFlow module to simulate 1-phase unsaturated flow, but I have encountered an issue. When I run the simulation using the PorousFlowUnsaturated Action versus not using the action, the results are significantly different.

Below are my input file and the simulation results. As you can see, there is a clear difference in pore pressure, and the results with the action seem more accurate.

Could someone please check if there is something wrong with my setup?

test1.i

# primary model units (m | s | Pa | kg --> N/m²)

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 10
    ymin = 0
    ymax = 6
    nx = 100
    ny = 60
  []
  [m1]
    input = mesh
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '10 2 0'
  []
  [injection_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'injection_area'
    normal = '-1 0 0'
    input = 'm1'
  []
  [outflow_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'outflow_area'
    normal = '1 0 0'
    input = 'injection_area'
  []
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
  PorousFlowDictator = dictator
[]

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

[Functions]
  [up_pressure_0]
    type = ParsedFunction
    expression = 10000*(6-y) #Pa
  []
  [down_pressure_0]
    type = ParsedFunction
    expression = 10000*(2-y) #Pa
  []
[]

[PorousFlowUnsaturated]
    porepressure = porepressure
    coupling_type = HydroMechanical
    gravity = '0 -10 0'
    fp = the_simple_fluid
    relative_permeability_exponent = 2
    relative_permeability_type = Corey
    residual_saturation = 0.05
    van_genuchten_alpha = 2E-3
    van_genuchten_m = 0.6
    stabilization = FULL
    flux_limiter_type = None
[]

[BCs]
  [archor_x]
    type = DirichletBC
    boundary = 'bottom'
    variable = disp_x
    value = 0
  []

  [archor_y]
    type = DirichletBC
    boundary = 'bottom'
    variable = disp_y
    value = 0
  []

  [Pressure]
    [Up_0]
      boundary = 'left'
      function = up_pressure_0
    []
    [Down_0]
      boundary = 'outflow_area'
      function = down_pressure_0
    []
  []

  [up_stream_pp_0]
    type = FunctionDirichletBC
    variable = porepressure
    function = up_pressure_0
    boundary = 'left'
  []
  [down_stream_pp_0]
    type = FunctionDirichletBC
    variable = porepressure
    function = down_pressure_0
    boundary = 'outflow_area'
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -10
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1E-3
    density0 = 1000.0
  []
[]

# ===== AuxVariable & AuxKernel: stress =====
[AuxVariables]
  [stress_maxprincipal]
    order = CONSTANT
    family = MONOMIAL

  []
  [stress_minprincipal]
    order = CONSTANT
    family = MONOMIAL

  []
  [saturation]
    family = MONOMIAL
    order = CONSTANT

  []
  [permeability]
    family = MONOMIAL
    order = CONSTANT

  []

  [porosity]
    order = CONSTANT
    family = MONOMIAL

  []

  [darcy_velocity_x]
    family = MONOMIAL
    order = CONSTANT

  []
  [darcy_velocity_y]
    family = MONOMIAL
    order = CONSTANT

  []
[]

[AuxKernels]
  [stress_maxprincipal]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_minprincipal
    scalar_type = MaxPrincipal
  []
  [stress_minprincipal]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_maxprincipal
    scalar_type = MinPrincipal
  []
  [saturation]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = saturation
  []
  [permeability]
    type = PorousFlowPropertyAux
    property = permeability
    column = 0
    row = 0
    variable = permeability
  []

  [porosity]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
  []

  [darcy_velocity_x]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_velocity_x
    component = x
    gravity = '0 -10 0'
  []
  [darcy_velocity_y]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_velocity_y
    component = y
    gravity = '0 -10 0'
  []
[]

# ===== Materials =====
[Materials]
  [strain]
    type = ComputeIncrementalStrain
  []

  [stress]
    type = ComputeFiniteStrainElasticStress
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 50e6
    poissons_ratio = 0.3
  []

  [density]
    type = GenericConstantMaterial
    prop_names = density
    prop_values = 2000
  []

  [permeability_val]
    type = PorousFlowPermeabilityConst
    permeability = '1e-10 0 0 0 1e-10 0 0 0 1e-10'
  []

  [porosity_val]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

# ===== Executioner =====
[Executioner]
  type = Transient
  automatic_scaling = true

  end_time = 10e3

  [TimeSteppers]
    [ConstDT1]
      type = FunctionDT
      function = 10
    []
  []

  solve_type = Newton

  petsc_options = '-snes_converged_reason'

  # best overall
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       mumps'

  nl_abs_tol = 1E-1
  nl_rel_tol = 1E-5
  nl_max_its = 200

  l_tol = 1E-10
  l_max_its = 200
[]

[Outputs]
  [out_vtk]
    type = VTK
    file_base = './test2_vtk/'
  []
  perf_graph = true
[]

[Debug]
  show_var_residual_norms = true
[]

test2.i

# primary model units (m | s | Pa | kg --> N/m²)

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 10
    ymin = 0
    ymax = 6
    nx = 100
    ny = 60
  []
  [m1]
    input = mesh
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '10 2 0'
  []
  [injection_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'injection_area'
    normal = '-1 0 0'
    input = 'm1'
  []
  [outflow_area]
    type = SideSetsAroundSubdomainGenerator
    block = 1
    new_boundary = 'outflow_area'
    normal = '1 0 0'
    input = 'injection_area'
  []
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
  use_displaced_mesh = false
  PorousFlowDictator = dictator
[]

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

[Functions]
  [up_pressure_0]
    type = ParsedFunction
    expression = 10000*(6-y) #Pa
  []
  [down_pressure_0]
    type = ParsedFunction
    expression = 10000*(2-y) #Pa
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_x disp_y'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc_rock]
    type = PorousFlowCapillaryPressureVG
    m = 0.60
    alpha = 2e-3 #Pa^-1
    sat_lr = 0.05
    # pc_max = 10e3
  []
[]

[BCs]
  [archor_x]
    type = DirichletBC
    boundary = 'bottom'
    variable = disp_x
    value = 0
  []

  [archor_y]
    type = DirichletBC
    boundary = 'bottom'
    variable = disp_y
    value = 0
  []

  [Pressure]
    [Up_0]
      boundary = 'left'
      function = up_pressure_0
    []
    [Down_0]
      boundary = 'outflow_area'
      function = down_pressure_0
    []
  []

  [up_stream_pp_0]
    type = FunctionDirichletBC
    variable = porepressure
    function = up_pressure_0
    boundary = 'left'
  []
  [down_stream_pp_0]
    type = FunctionDirichletBC
    variable = porepressure
    function = down_pressure_0
    boundary = 'outflow_area'
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -10
  []

  [stress_x]
    type = StressDivergenceTensors
    component = 0
    variable = disp_x
  []
  [stress_y]
    type = StressDivergenceTensors
    component = 1
    variable = disp_y
  []
  [effective_stress_x]
    type = PorousFlowEffectiveStressCoupling
    variable = 'disp_x'
    component = 0
  []

  [effective_stress_y]
    type = PorousFlowEffectiveStressCoupling
    variable = 'disp_y'
    component = 1
  []

  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = 'porepressure'
  []

  [flux]
    type = PorousFlowAdvectiveFlux
    use_displaced_mesh = false
    variable = porepressure
    gravity = '0 -10 0'
    fluid_component = 0
  []
  [poro_vol_exp]
    type = PorousFlowMassVolumetricExpansion
    variable = 'porepressure'
    fluid_component = 0
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 1E-3
    density0 = 1000.0
  []
[]

# ===== AuxVariable & AuxKernel: stress =====
[AuxVariables]
  [stress_maxprincipal]
    order = CONSTANT
    family = MONOMIAL

  []
  [stress_minprincipal]
    order = CONSTANT
    family = MONOMIAL

  []
  [saturation]
    family = MONOMIAL
    order = CONSTANT

  []
  [permeability]
    family = MONOMIAL
    order = CONSTANT

  []

  [porosity]
    order = CONSTANT
    family = MONOMIAL

  []

  [darcy_velocity_x]
    family = MONOMIAL
    order = CONSTANT

  []
  [darcy_velocity_y]
    family = MONOMIAL
    order = CONSTANT

  []
[]

[AuxKernels]
  [stress_maxprincipal]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_minprincipal
    scalar_type = MaxPrincipal
  []
  [stress_minprincipal]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_maxprincipal
    scalar_type = MinPrincipal
  []
  [saturation]
    type = MaterialStdVectorAux
    property = PorousFlow_saturation_qp
    index = 0
    variable = saturation
  []
  [permeability]
    type = PorousFlowPropertyAux
    property = permeability
    column = 0
    row = 0
    variable = permeability
  []

  [porosity]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
  []

  [darcy_velocity_x]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_velocity_x
    component = x
    gravity = '0 -10 0'
  []
  [darcy_velocity_y]
    type = PorousFlowDarcyVelocityComponent
    variable = darcy_velocity_y
    component = y
    gravity = '0 -10 0'
  []
[]

# ===== Materials (linear-elastic to keep it simple) =====
[Materials]
  [strain]
    type = ComputeIncrementalStrain
  []

  [stress]
    type = ComputeFiniteStrainElasticStress
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 50e6
    poissons_ratio = 0.3
  []

  [density]
    type = GenericConstantMaterial
    prop_names = density
    prop_values = 2000
  []

  [permeability_val]
    type = PorousFlowPermeabilityConst
    permeability = '1e-10 0 0 0 1e-10 0 0 0 1e-10'
  []

  [porosity_val]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []

  [temperature]
    type = PorousFlowTemperature
  []
  [eff_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [vol_strain]
    type = PorousFlowVolumetricStrain
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = the_simple_fluid
    phase = 0
  []

  [ppss_rock]
    type = PorousFlow1PhaseP
    porepressure = porepressure
    capillary_pressure = pc_rock
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.05
    sum_s_res = 0.05
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

# ===== Executioner =====
[Executioner]
  type = Transient
  automatic_scaling = true

  end_time = 10e3

  [TimeSteppers]
    [ConstDT1]
      type = FunctionDT
      function = 10
    []
  []

  solve_type = Newton

  petsc_options = '-snes_converged_reason'

  # best overall
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       mumps'

  nl_abs_tol = 1E-1
  nl_rel_tol = 1E-5
  nl_max_its = 200

  l_tol = 1E-10
  l_max_its = 200
[]

[Outputs]
  [out_vtk]
    type = VTK
    file_base = './test2_2_vtk/'
  []
  perf_graph = true
[]

[Debug]
  show_var_residual_norms = true
[]

[GiudGiud]

Hello

Use Problem/type=DumpObjectsProblem on both inputs to see what objects the Action is creation and how the two do not match

Guillaume

[PEI0214]

Thank you for your help. I have found the difference. However, I’ve encountered another issue: setting pc_max doesn’t seem to limit the maximum negative pore pressure. For example, when I set the maximum matric suction to 10 kPa and 100 kPa, the results were identical. Moreover, the negative pore pressure greatly exceeds the specified pc_max. Why is this happening?

[GiudGiud]

I m not familiar with this. How is pc_max enforced? Using the Bounds system? If not you could try that