Point Source and PolyLine Dirac production convergence issues

[Stevenovia]

See my .i file below

The exact injection scenario for the simulation runs very well but when the Dirac source has their sign flipped for production, then it becomes an issue as it never reaches the simulation end. It begins to have convergence issues before 100 seconds.
The .i file is attached below.

PORO_1MPa_Initial_Cond.txt

Also see the details below

[Mesh]
  type = FileMesh
  file = GW_30_07.exo
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
  displacements = 'disp_x disp_y disp_z'
[]

[Variables]
  [porepressure]
    initial_condition = 1e6    # 1 MPa
    scaling = 1e-6
  []
  [disp_x]
    scaling = 1e-10
  []
  [disp_y]
    scaling = 1e-10
  []
  [disp_z]
    scaling = 1e-10
  []
[]
[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
[]

# --- Essential Kernels only ---
[Kernels]
  # Flow: storage + Darcy flux + coupling to volumetric strain
  [mass_time]
    type = PorousFlowMassTimeDerivative
    variable = porepressure
  []
  [darcy]
    type = PorousFlowFullySaturatedDarcyFlow
    variable = porepressure
  []
  [vol_strain_cpl]
    type = PorousFlowMassVolumetricExpansion
    variable = porepressure
  []

  # Mechanics: divergence of stress in each direction
  [mech_x]
    type = StressDivergenceTensors
    variable = disp_x
    component = 0
    use_displaced_mesh = false
  []
  [mech_y]
    type = StressDivergenceTensors
    variable = disp_y
    component = 1
    use_displaced_mesh = false
  []
  [mech_z]
    type = StressDivergenceTensors
    variable = disp_z
    component = 2
    use_displaced_mesh = false
  []

  # Poro-mechanical effective stress coupling (α=1 here; adjust if needed)
  [poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    component = 0
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    component = 1
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_z
    component = 2
  []
[]

# --- Mechanical BCs only (no pressure BCs = no-flow boundaries) ---
[BCs]
  # X faces fixed in X
  [xneg_fixed_x]
    type = DirichletBC
    variable = disp_x
    boundary = boundary_xneg
    value = 0
  []
  [xpos_fixed_x]
    type = DirichletBC
    variable = disp_x
    boundary = boundary_xpos
    value = 0
  []

  # Y faces fixed in Y
  [yneg_fixed_y]
    type = DirichletBC
    variable = disp_y
    boundary = boundary_yneg
    value = 0
  []
  [ypos_fixed_y]
    type = DirichletBC
    variable = disp_y
    boundary = boundary_ypos
    value = 0
  []

  # Bottom fixed in Z; top free (stress-free)
  [zneg_fixed_z]
    type = DirichletBC
    variable = disp_z
    boundary = boundary_zneg
    value = 0
  []
  # Porepressure fixed to 0 on all side boundaries
  [p_all_zero]
    type = DirichletBC
    variable = porepressure
    boundary = 'boundary_xneg boundary_xpos boundary_yneg boundary_ypos boundary_zneg'
    value = 1e6
  []
[]

# --- Fluid properties (simple water) ---
[FluidProperties]
  [water]
    type = SimpleFluidProperties
    bulk_modulus = 2.2e9
    density0 = 1000
    viscosity = 1e-3
  []
[]

# --- Minimal materials per layer (elasticity, small strain, stress, perm, φ const, Biot, effective pressure) ---
[Materials]
  # Core PorousFlow materials
  [pp_sat]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  
  []
  # --- Added: provide temperature (QP & nodal) on all rock blocks ---
  [T_qp]
    type = PorousFlowTemperature
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  [T_nodal]
    type = PorousFlowTemperature
    at_nodes = true
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  # --- Added: provide mass fraction (QP & nodal) for single-component ---
  [mass_frac]
    type = PorousFlowMassFraction
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  [mass_frac_nodal]
    type = PorousFlowMassFraction
    at_nodes = true
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  [fluid]
    type = PorousFlowSingleComponentFluid
    fp = water
    phase = 0
  []
  [p_eff]
    type = PorousFlowEffectiveFluidPressure
  []
  [vstrain]
    type = PorousFlowVolumetricStrain
  []
  
  # BASIN_FILL Layer Properties
  [elasticity_basin_fill]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 2.198e10  # Pa
    poissons_ratio = 0.34
    block = Basin_Fill
  []
  [strain_basin_fill]
    type = ComputeSmallStrain
    block = Basin_Fill
  []
  [stress_basin_fill]
    type = ComputeLinearElasticStress
    block = Basin_Fill
  []
  [porosity_basin_fill]
    type = PorousFlowPorosityConst
    porosity = 0.12  # 12%
    block = Basin_Fill
  []
  [permeability_basin_fill]
    type = PorousFlowPermeabilityConst
    permeability = '1.7e-14 0 0 0 1.7e-14 0 0 0 1.7e-14'  # m^2
    block = Basin_Fill
  []
  [biot_basin_fill]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 1.0
    solid_bulk_compliance = 4.55e-11  # 1/youngs_modulus approximately
    block = Basin_Fill
  []
  
  # CLAY Layer Properties  
  [elasticity_clay]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.202e10  # Pa
    poissons_ratio = 0.333
    block = Clay
  []
  [strain_clay]
    type = ComputeSmallStrain
    block = Clay
  []
  [stress_clay]
    type = ComputeLinearElasticStress
    block = Clay
  []
  [porosity_clay]
    type = PorousFlowPorosityConst
    porosity = 0.07  # 7%
    block = Clay
  []
  [permeability_clay]
    type = PorousFlowPermeabilityConst
    permeability = '1.0e-16 0 0 0 1.0e-16 0 0 0 1.0e-16'  # m^2
    block = Clay
  []
  [biot_clay]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 1.0
    solid_bulk_compliance = 8.32e-11  # 1/youngs_modulus approximately
    block = Clay
  []
  
  # VOLCANICS Layer Properties
  [elasticity_volcanics]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 4.916e10  # Pa
    poissons_ratio = 0.256
    block = Volcanics
  []
  [strain_volcanics]
    type = ComputeSmallStrain
    block = Volcanics
  []
  [stress_volcanics]
    type = ComputeLinearElasticStress
    block = Volcanics
  []
  [porosity_volcanics]
    type = PorousFlowPorosityConst
    porosity = 0.07  # 7%
    block = Volcanics
  []
  [permeability_volcanics]
    type = PorousFlowPermeabilityConst
    permeability = '1.0e-17 0 0 0 1.0e-17 0 0 0 1.0e-17'  # m^2
    block = Volcanics
  []
  [biot_volcanics]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 1.0
    solid_bulk_compliance = 2.03e-11  # 1/youngs_modulus approximately
    block = Volcanics
  []
  
  # BASEMENT Layer Properties
  [elasticity_basement]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 7.028e10  # Pa
    poissons_ratio = 0.234
    block = Basement
  []
  [strain_basement]
    type = ComputeSmallStrain
    block = Basement
  []
  [stress_basement]
    type = ComputeLinearElasticStress
    block = Basement
  []
  [porosity_basement]
    type = PorousFlowPorosityConst
    porosity = 0.01  # 1%
    block = Basement
  []
  [permeability_basement]
    type = PorousFlowPermeabilityConst
    permeability = '1.0e-18 0 0 0 1.0e-18 0 0 0 1.0e-18'  # m^2
    block = Basement
  []
  [biot_basement]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 1.0
    solid_bulk_compliance = 1.42e-11  # 1/youngs_modulus approximately
    block = Basement
  []
[]


[AuxVariables]
  # element-wise normal strains
  [eps_xx]
    family = MONOMIAL
    order = CONSTANT
  []
  [eps_yy]
    family = MONOMIAL
    order = CONSTANT
  []
  [eps_zz]
    family = MONOMIAL
    order = CONSTANT
  []

  # element-wise principal stresses
  [sigma_max]
    family = MONOMIAL
    order = CONSTANT
  []
  [sigma_mid]
    family = MONOMIAL
    order = CONSTANT
  []
  [sigma_min]
    family = MONOMIAL
    order = CONSTANT
  []
[]


[AuxKernels]
  [eps_xx_aux]
    type = RankTwoAux
    variable = eps_xx
    rank_two_tensor = mechanical_strain
    index_i = 0
    index_j = 0
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  [eps_yy_aux]
    type = RankTwoAux
    variable = eps_yy
    rank_two_tensor = mechanical_strain
    index_i = 1
    index_j = 1
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  [eps_zz_aux]
    type = RankTwoAux
    variable = eps_zz
    rank_two_tensor = mechanical_strain
    index_i = 2
    index_j = 2
    block = 'Basin_Fill Clay Volcanics Basement'
  []

  # principal stresses (note the allowed enums)
  [sig_max_aux]
    type = RankTwoScalarAux
    variable = sigma_max
    rank_two_tensor = stress
    scalar_type = MaxPrincipal
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  [sig_mid_aux]
    type = RankTwoScalarAux
    variable = sigma_mid
    rank_two_tensor = stress
    scalar_type = MidPrincipal
    block = 'Basin_Fill Clay Volcanics Basement'
  []
  [sig_min_aux]
    type = RankTwoScalarAux
    variable = sigma_min
    rank_two_tensor = stress
    scalar_type = MinPrincipal
    block = 'Basin_Fill Clay Volcanics Basement'
  []
[]


[Functions]
  [well_ramp]
    type = PiecewiseLinear
    # time [s] :    0     1 hr   3 hr
    x = '0  3600  10800'
    y = '0   0.5     1'
  []
[]

[DiracKernels]
  [inj_force1]
    type = PorousFlowSquarePulsePointSource
    variable = porepressure
    point = '334.200 4263.000 1.466'
    mass_flux = 8.2          # + injects, - produces
    start_time = 0
    end_time = 3600          # match your Executioner end_time
  []
  [inj_force2]
    type = PorousFlowSquarePulsePointSource
    variable = porepressure
    point = '334.200 4263.000 1.466'
    mass_flux = 41.00          # + injects, - produces
    start_time = 3600
    end_time = 8.64e6          # match your Executioner end_time
  []
  [inj_force3]
    type = PorousFlowSquarePulsePointSource
    variable = porepressure
    point = '333.600 4263.000 1.449'
    mass_flux = 3.154
    start_time = 0
    end_time = 3600
  []
  [inj_force4]
    type = PorousFlowSquarePulsePointSource
    variable = porepressure
    point = '333.600 4263.000 1.449'
    mass_flux = 15.77
    start_time = 3600
    end_time = 8.64e6
  []
  [inj_force5]
    type = PorousFlowSquarePulsePointSource
    variable= porepressure
    point = '332.700 4263.000 0.987'
    mass_flux = 12.618
    start_time = 0
    end_time = 3600
  []
  [inj_force6]
    type = PorousFlowSquarePulsePointSource
    variable = porepressure
    point = '332.700 4263.000 0.987'
    mass_flux = 63.09
    start_time = 3600
    end_time = 8.64e6
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_hypre_type'
    petsc_options_value = 'hypre boomeramg'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  line_search = bt

  start_time = 0
  end_time = 8.64e6

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 10                 # start at 10 s
    optimal_iterations = 8  # target nonlinear its
    growth_factor = 1.3
    cutback_factor = 0.5
    
  []

  # tolerances: moderate nonlinear, not-insanely-tight linear
  nl_rel_tol = 1e-3
  nl_abs_tol = 1e-8
  nl_max_its = 30
  l_tol = 1e-6
  l_max_its = 200
[]


[Outputs]
  execute_on = 'initial timestep_end'
  exodus = true
  csv = true
[]

[GiudGiud]

Hello

Please see these instructions for figuring out convergence issues
https://mooseframework.inl.gov/application_usage/failed_solves.html

nl_rel_tol = 1e-3

This seems too lose. Especially with real low scaling factors

scaling = 1e-6/1e-10

So maybe try removing these settings and using automatic scaling?

[Stevenovia]

I have tried using them but no improvements. The scaling was changed to automatic scaling. The problem is that for injection, it works very well, but for production it just has convergence problems right from the first time step.

[GiudGiud]

can you ramp up production from a low value or does it always not converge?
What did the SVD monitor say on a low-DOF case?

[Stevenovia]

I also ramped up production from a low value but the result was the same. It always does not converge.

I later changed the setup a bit and used a polyline that also had production ramped up slowly. See the Diracs below for that case. This one stalled at about 115000 seconds. The pressure dropped at the bottomhole (at 3342) to a max of about -9 x 10^9 Pascals when a stall was experienced. By stall I mean that the simulation started having convergence issues and the dt was then slowly cutback until it just could not go any further.

# --- time ramp for gentle turn-on of injection ---
[Functions]
  [well_ramp]
    type = PiecewiseLinear
    x = '0   3600   10800'   # s
    y = '0   0.2    1.0'
  []
[]
[AuxVariables]
  [ramp]
    family = MONOMIAL
    order = CONSTANT
  []
[]
[AuxKernels]
  [ramp_aux]
    type = FunctionAux
    variable = ramp
    function = well_ramp
    execute_on = 'initial timestep_begin'
  []
[]

# --- INJECTION via PolyLineSink (constant total rate * ramp) ---
[DiracKernels]
  [prod_3327]
    type = PorousFlowPolyLineSink
    variable = porepressure
    SumQuantityUO = inj_mass_3327
    point_file = well_3327.bh
    function_of = pressure
    fluid_phase = 0
    # floor ……… mid ………… near-full ………… full (Pa)
    p_or_t_vals = '6.0e5   8.5e5          9.6e5          1.0e6'
    # fluxes [kg/(m·s)] at those pressures (keep last = your original)
    fluxes      = '8       20.5           32              41'
    use_mobility = false
    multiplying_var = ramp
  []
  [prod_3336]
    type = PorousFlowPolyLineSink
    variable = porepressure
    SumQuantityUO = inj_mass_3336
    point_file = well_3336.bh
    function_of = pressure
    fluid_phase = 0
    p_or_t_vals = '6.0e5   8.5e5          9.6e5          1.0e6'
    fluxes      = '3       7.5            12              15'
    use_mobility = false
    multiplying_var = ramp
  []
  [prod_3342]
    type = PorousFlowPolyLineSink
    variable = porepressure
    SumQuantityUO = inj_mass_3342
    point_file = well_3342.bh
    function_of = pressure
    fluid_phase = 0
    p_or_t_vals = '6.0e5   8.5e5          9.6e5          1.0e6'
    fluxes      = '12      30.5           48              61'
    use_mobility = false
    multiplying_var = ramp
  []
[]

[GiudGiud]

can you please paste the simulation log? To see how it was failing.
@rpodgorney or @cpgr could have an idea why this would fail?

[rpodgorney]

Place a post-processor at the injection locations, send along a plot of pressure, density, porosity, permeability, etc. I have some thoughts but would like to see these data before I comment