Residual saturation appears to be ignored in PorousFlow

[sinaomrani96]

Dear all,

I hope you are well.

I am trying to simulate a two-phase, three-component system consisting of CO₂, brine, and salt. The simulation initially runs as expected, but at a certain point it encounters convergence problems.

To investigate the issue, I simplified my model using the existing PorousFlow examples. However, I still could not run the simulation long enough for the domain to become fully saturated with CO₂. I am currently trying to identify the cause of the convergence failure. The error that I get if I run the simulation long enough is:

21 Nonlinear |R| = 1.306216e-06
      0 Linear |R| = 1.306216e-06
      1 Linear |R| = 9.013093e-21
    |residual|_2 of individual variables:
                      pgas:  7.10773e-07
                      z:     9.39817e-07
                      xnacl: 3.60454e-07
22 Nonlinear |R| = 1.232226e-06
    Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0
                   PC failed due to FACTOR_NUMERIC_ZEROPIVOT

Temperature -nan is out of range in water: inRegion()
To recover, the solution will fail and then be re-attempted with a reduced time step.

  Nonlinear solve did not converge due to DIVERGED_FUNCTION_NANORINF iterations 22
  Finished Solving                                                                       [ 30.17 s] [  193 MB]
 Solve Did NOT Converge!
Aborting as solve did not converge

I have read PorousFlow documents regarding issues with near Seff 0 and 1. But still think this is not normal. While examining the results, I noticed two potentially unexpected behaviours:

1. The gas saturation appears to exceed the value allowed by the residual water saturation. In other words, the water saturation falls below the specified residual water saturation.
2. Near the injection point, the CO2 saturation begins to decrease despite continuous CO2 injection.

I tested several capillary-pressure and relative-permeability formulations and parameter values, but these behaviours persisted. I am therefore unsure whether this is related to the kernels, BCs, the nonlinear/linear solver settings, or a possible issue in my input file.

I have included the input file below, along with figures showing the saturation behaviour and the results near the injection boundary.

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 40
    ny = 10
    xmin = 0
    xmax = 4
    ymin = 0
    ymax = 4
  []
  coord_type = RZ
  rz_coord_axis = Y
[]


[GlobalParams]
  PorousFlowDictator = dictator
  temperature = 50
  gravity = '0 0 0'
[]

[Variables]
  [pgas]
    initial_condition = 20e6
  []
  [z]
    initial_condition = 0
    # scaling = 1e3
  []
  [xnacl]
    initial_condition = 0.1
  []
[]

[DiracKernels]
  [source]
    type = PorousFlowSquarePulsePointSource
    variable = z
    point = '0 2.2 0'
    mass_flux = 0.0001
  []
[]


[BCs]
  [right_gas_component2]
    type = PorousFlowOutflowBC
    boundary = right
    variable = z
    mass_fraction_component = 1
  []
  [pgas_l]
    type = DirichletBC
    boundary = right
    variable = pgas
    value = 20e6
  []
  [salt_r]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = xnacl
    multipliers = '0 1e9'
    pt_vals = '0 1e9'
    fluid_phase = 0
    flux_function = 1
    use_mobility = true
    use_relperm = true
    mass_fraction_component = 2
  []
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'initial timestep_end'
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    variable = z
    fluid_component = 1
  []
  [mass2]
    type = PorousFlowMassTimeDerivative
    variable = xnacl
    fluid_component = 2
  []
  [flux2]
    type = PorousFlowAdvectiveFlux
    variable = xnacl
    fluid_component = 2
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z xnacl'
    number_fluid_phases = 2
    number_fluid_components = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
    sat_lr = 0.2
  []
  # [pc]
  #   type = PorousFlowCapillaryPressureVG
  #   m = 0.9
  #   alpha = 1e-2
  #   pc_max = 1e8
  #   s_scale = 0.5
  # []
  # [pc]
  #   type = PorousFlowCapillaryPressureBC
  #   lambda = 0.7
  #   pe = 1e4
  # []
  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
  [brine]
    type = BrineFluidPropertiesBeta
    water_fp = water
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  # [relperm0]
  #   type = PorousFlowRelativePermeabilityCorey
  #   n = 2
  #   phase = 0
  #   s_res = 0.2
  #   sum_s_res = 0.3
  # []
  # [relperm1]
  #   type = PorousFlowRelativePermeabilityCorey
  #   n = 3
  #   phase = 1
  #   s_res = 0.1
  #   sum_s_res = 0.3
  # []
  [relperm0]
    type = PorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = 0.2
    sum_s_res = 0.3
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityBC
    phase = 1
    lambda = 2
    nw_phase = true
    s_res = 0.1
    sum_s_res = 0.3
  []

[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1e10
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    growth_factor = 1.5
  []
  # controls for linear iterations
  l_max_its = 40
  l_abs_tol = 1e-20
  # controls for nonlinear iterations
  nl_max_its = 40
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8
[]

[Preconditioning]
  active = mumps_is_best_for_parallel_jobs
  [mumps_is_best_for_parallel_jobs]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_factor_shift_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu   NONZERO    mumps'
  []
  [lu_mumps]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -pc_factor_shift_type -pc_factor_shift_amount'
    petsc_options_value = 'lu       mumps                         NONZERO               1e-10'
  []
  [basic]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2             '
  []
  [preferred]
        type = SMP
        full = true
        petsc_options_iname = '-pc_type -sub_pc_factor_shift_type -pc_factor_mat_solver_package -snes_type -pc_factor_shift_type -snes_linesearch_type'
        petsc_options_value = 'lu NONZERO mumps vinewtonssls nonzero l2'
  []
  [preferred1]
        type = SMP
        full = true
        petsc_options_iname = '-pc_type -sub_pc_type -pc_factor_mat_solver_package -snes_type -pc_factor_shift_type -snes_linesearch_type'
        petsc_options_value = 'bjacobi lu mumps vinewtonssls NONZERO l2'
  []
  [smp1]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
    petsc_options_value = 'gmres bjacobi lu NONZERO'
  []
  [smp]
    # this block is part of what is being tested, see "tests" file
    type = SMP
    full = true
  []
[]

[Outputs]
  print_linear_residuals = true
  perf_graph = true
  exodus = true
  [csvout]
    type = CSV
    execute_on = timestep_end
    execute_vector_postprocessors_on = final
  []
[]

[Debug]
  show_var_residual_norms = true
     # show_top_residuals = 2
     # show_execution_order = 'INITIAL TIMESTEP_BEGIN'
[]

As you see below sat_gas goes to 0.97 while s_res=0.2. I would greatly appreciate any comments or suggestions on what may be causing these convergence and saturation issues.

Many thanks for your help.

[cpgr]

This CO2-brine EOS allows the water phase to evaporate into the gas phase, as well as the CO2 to dissolve into the water. So you find that even though liquid rel perm is 0 for Sl < 0.2, and it can't advect away, it is evaporated into the gas phase, where it can flow out.

This drying out behaviour is definitely a feature!

On the convergence - this is just hard. I've had success with line_search = basic when it's been really slow convergence. You could try that. Feel free to email me to discuss this work too.