Multiphase Porous Flow (pliquid/h) does not converge with steam

[aikubo]

Hello all (@cpgr especially)

I am working on modeling the hydrothermal system of a dike intruded into granite and struggling with nonlinear convergence. I've prepared a simple example to demonstrate but the goal is to eventually

1. benchmark with USGS Hydrotherm Code
2. couple this to a Navier Stokes flow model multiapp of flow of magma in the dike
3. conjugate heat transfer with magma in the dike cooling

The problem is as follows:
A dike intrudes into granite host rock at temperature of 500 C. It is a constant temperature BC for 2 (6.3e7 s) then shuts off and cools (Neuman BC, value = 0).
The boundary at the dike is impermeable to flow (NeumanBC for porepressure).
I've expanded the domain to include what I consider the "far field" so the right most boundaries
should be uneffected by the dike heating (4000m).
The simulation stops at 200 m depth and goes down to 2000 m depth. There is a hydrostatic pressure gradient and geothermal temperature gradient. Water is the only fluid phase.

So far to deal with convergence issues I've tried the following

1. Using PorousFlowPiecewiseLinearSink BCs instead of DirichletBC - This helped it run longer
2. Using SimpleFluidProperties instead of Water97FluidProperties - This helps the simulation runs to completion and works fine but I'd like to use water properties and especially want to look at multiphase effects and benchmarking
3. Changing permeablitiy - Decreasing permeability to 1e-15 causes it to fail to converge after only a few timesteps
4. Trying different BCs such as `PorousFlowOutflowBC' didn't help.
5. Adding dampers to pressure and temperature. This helps a bit but not in the long run since it will eventually cause convergence issues becuase I keep getting unphysical pressure.
6. Ramping up temperature (enthalpy) slowly over time - it stops converging when the temperature reaches a critical point and phase chanage occurs ~5000-10000 seconds
7. constant capillary pressure with pc = 0
8. If i lower the enthalpy BC to a level with no phase change it converges fine and runs to completion
9. reducing nl_abs_tol to as low as 1e-4, still convergence issues associated with steam

I believe this is a boundary condition issue because when I get problems with pressure or temperature out of bounds (dampers off) it usually occurs on the top right corner or left corner.
Is this boundary conditions, initial conditions (try an IC with some steam initially?), or maybe capillary pressure?

Things I've observed that seem strange to me :

1. despite scaling the residual for enthalpy "h" is often much larger than for pressure and increases faster, Does it need more BCS?
2. simulation reaches a point where there is no change in residual despite iterations for example:

[DBG][0] Max 1 residuals
[DBG][0] -2191180.71139573 'h' in subdomain(s) {''} at node 441: (x,y,z)=(     100,     -380,        0)
 6 Nonlinear |R| = 2.597688e+06
      0 Linear |R| = 2.597688e+06
      1 Linear |R| = 1.601581e-09
 Damping factor: 0
    |residual|_2 of individual variables:
                    pliquid: 25.6896
                    h:       2.59769e+06
[DBG][0] Max 1 residuals
[DBG][0] -2191180.71139573 'h' in subdomain(s) {''} at node 441: (x,y,z)=(     100,     -380,        0)
 7 Nonlinear |R| = 2.597688e+06
      0 Linear |R| = 2.597688e+06
      1 Linear |R| = 1.601581e-09
 Damping factor: 0
    |residual|_2 of individual variables:
                    pliquid: 25.6896
                    h:       2.59769e+06
[DBG][0] Max 1 residuals
[DBG][0] -2191180.71139573 'h' in subdomain(s) {''} at node 441: (x,y,z)=(     100,     -380,        0)
 8 Nonlinear |R| = 2.597688e+06
      0 Linear |R| = 2.597688e+06
      1 Linear |R| = 1.601581e-09
 Damping factor: 0
    |residual|_2 of individual variables:
                    pliquid: 25.6896
                    h:       2.59769e+06
[DBG][0] Max 1 residuals
[DBG][0] -2191180.71139573 'h' in subdomain(s) {''} at node 441: (x,y,z)=(     100,     -380,        0)
 9 Nonlinear |R| = 2.597688e+06
      0 Linear |R| = 2.597688e+06
      1 Linear |R| = 1.601581e-09

Here is a simplified version of the input file:

# can you have an adaptive BC?
# perfectly matched layers 
# check the pressures at the BC and the flux through the boundary 
# doesn't have to be atmospheric

# nusselt number 
# scaling 
# heat flux at surface - can caluculate it from space 
# what is the thermal annomaly from space?

#trying nonzero cap pressure 
# might have issues because cap pressure is infinite when sat water is low 
# BC pressure is what is used in the watervapor.i model


[Mesh]
    [gen]
      type = GeneratedMeshGenerator
      dim = 2
      nx = 20
      ny = 20
      xmin = 0
      xmax = 4000
      ymax = -200
      ymin = -2000
    []
    [dike]
      type = ParsedGenerateSideset
      input = gen
      combinatorial_geometry = 'y <= -400 & x = 0'
      new_sideset_name = dike
    []

  []
  
  [Dampers]
    [./limit]
      type = BoundingValueNodalDamper
      variable = pliquid
      max_value = 1e8
      min_value = 1e1
    [../]
    [./limit2]
        type = BoundingValueNodalDamper
        variable = h
        max_value = 1e6
        min_value = 1e2

    [../]
  []

  [GlobalParams]
    PorousFlowDictator = dictator
    gravity = '0 -9.81 0'
  []

  [UserObjects]
    [dictator]
      type = PorousFlowDictator
      porous_flow_vars = 'pliquid h'
      number_fluid_phases = 2
      number_fluid_components = 1
    []
    [pc] #from porousflow/test/tests/fluidstate/watervapor.i
      type = PorousFlowCapillaryPressureBC
      pe = 1e5
      lambda = 2
      pc_max = 1e6
    []
    [fs]
      type = PorousFlowWaterVapor
      water_fp = water97 # you can't use tabulated fluids here!
      capillary_pressure = pc
    []
  []
  
  
  [AuxVariables]
    [temperature]
      family = MONOMIAL
      order = CONSTANT
    []
    [water_darcy_vel_x]
      family = MONOMIAL
      order = CONSTANT
    []
    [water_darcy_vel_y]
      family = MONOMIAL
      order = CONSTANT
    []
    [gas_darcy_vel_x]
      family = MONOMIAL
      order = CONSTANT
    []
    [gas_darcy_vel_y]
      family = MONOMIAL
      order = CONSTANT
    []
    [pgas]
      family = MONOMIAL
      order = CONSTANT
    []
    [gas_sat]
        family = MONOMIAL
        order = CONSTANT  
    []
    [hydrostat]
        family = MONOMIAL
        order = CONSTANT
    []
    [geotherm]
      order = FIRST
      family = LAGRANGE
    []
[]
  
  
  [AuxKernels]
    [temperature]
      type = PorousFlowPropertyAux
      variable = temperature
      property = temperature
      execute_on = 'initial timestep_end'  
    []
    [darcy_vel_x_kernel]
      type = PorousFlowDarcyVelocityComponent
      component = x
      variable = water_darcy_vel_x
      fluid_phase = 0                            
      execute_on = 'initial timestep_end'
       
    []
    [darcy_vel_y_kernel]
      type = PorousFlowDarcyVelocityComponent
      component = y
      variable = water_darcy_vel_y
      fluid_phase = 0                             
      execute_on = 'initial timestep_end'
    []
    [darcy_vel_x_kernel_gas]
      type = PorousFlowDarcyVelocityComponent
      component = x
      variable =gas_darcy_vel_x
      fluid_phase = 1                            
      execute_on = 'initial timestep_end'
       
    []
    [darcy_vel_y_kernel_gas]
      type = PorousFlowDarcyVelocityComponent
      component = y
      variable = gas_darcy_vel_y
      fluid_phase = 1                             
      execute_on = 'initial timestep_end'
       
    []
    [pressure_gas]
      type = PorousFlowPropertyAux
      variable = pgas
      property = pressure
      phase = 1
      execute_on = 'initial timestep_end'
    []
    [gas_sat]
        type = PorousFlowPropertyAux
        variable = gas_sat
        property = saturation
        phase = 1
        execute_on = 'initial timestep_end'
    []
    [hydrostat]
        type = FunctionAux
        function = ppfunc
        variable = hydrostat
    []
  []
  
  [Variables]
    [pliquid]
      order = FIRST
      family = LAGRANGE
      scaling = 1e1
    []
    [h]
      order = FIRST
      family = LAGRANGE
    []
  []
  
  [Functions]
    [dike_temp]
        type = ParsedFunction
        expression = '500000*(1-exp(-t/5000))+40000'    
      []
    [ppfunc]
      type = ParsedFunction
      expression = 1.0135e5-(y)*9.81*1000 #hydrostatic gradientose   + atmospheric pressure in Pa
    []
    [tfunc]
      type = ParsedFunction
      expression = 285+(-y)*10/1000 # geothermal 10 C per kilometer in kelvin
    []
  []
  

  
  [ICs]
    [t_ic]
      type = FunctionIC
      variable = geotherm
      function = tfunc
    []
    [ppic]
      type = FunctionIC # pressure is hydrostatic
      variable = pliquid
      function = ppfunc
    []
    [hic]
      type = PorousFlowFluidPropertyIC
      variable = h
      temperature = geotherm
      property = enthalpy
      temperature_unit = Kelvin
      porepressure = pliquid
      fp = water97
    []

  []
  
  [BCs]
    [t_dike_dirichlet]
        type = FunctionDirichletBC
        variable = h
        function = dike_temp
        boundary = 'dike'
    []
    [t_dike_neumann]
        type = NeumannBC
        variable = h
        boundary = 'dike'
        value = 0
    []

    [pp_like_dirichlet]
        type = PorousFlowPiecewiseLinearSink
        variable = pliquid
        boundary = 'top bottom right'
        pt_vals = '1e-9 1e9'
        multipliers = '1e-9 1e9'
        PT_shift = hydrostat
        flux_function = 1e-5
        use_mobility = true 
        use_relperm = true
        fluid_phase = 0
    []
    [ptop_gas]
        type = PorousFlowPiecewiseLinearSink
        # allow fluid to flow out or in of the top boundary
        # based on pliquid - Pe
        variable = pliquid
        boundary = 'top bottom right'
        pt_vals = '1e-9 1e9'
        multipliers = '1e-9 1e9'
        PT_shift = hydrostat
        flux_function = 1e-5
        fluid_phase = 1
        use_mobility = true
        use_relperm = true
    []
  []
  
  [Kernels]
    [mass]
      type = PorousFlowMassTimeDerivative
      variable = pliquid
      
    []
    [massflux]
      type = PorousFlowAdvectiveFlux
      variable = pliquid
       
    []
    [heat]
      type = PorousFlowEnergyTimeDerivative
      variable = h
       
    []
    [heatflux]
      type = PorousFlowHeatAdvection
      variable = h
       
    []
    [heatcond]
      type = PorousFlowHeatConduction
      variable = h 
    []
  []
  
 
  [FluidProperties]
    [water97]
      type = Water97FluidProperties    # IAPWS-IF97
    []
    [water_tab]
      type = TabulatedBicubicFluidProperties
      fp = water97
      temperature_min=274
      temperature_max=1000
      pressure_min=1
      pressure_max=1e9
      error_on_out_of_bounds = false
      save_file = true
    []
  []
  
  [Materials]
      [watervapor]
        type = PorousFlowFluidStateSingleComponent
        porepressure = pliquid
        enthalpy = h
        capillary_pressure = pc
        fluid_state = fs
        temperature_unit = Kelvin
      []

      [porosity_wallrock]
        type = PorousFlowPorosityConst
        porosity = 0.1
      []
      [permeability]
        type = PorousFlowPermeabilityConst
        permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-15'
      []
      [relperm_water] # from watervapor.i
        type = PorousFlowRelativePermeabilityCorey
        n = 2
        phase = 0
      []
      [relperm_gas]  # from watervapor.i
        type = PorousFlowRelativePermeabilityCorey
        n = 3
        phase = 1
      []
      [internal_energy] 
        type = PorousFlowMatrixInternalEnergy
        density = 2400 # kg/m^3
        specific_heat_capacity = 790 # J/kg/K
      []
      [rock_thermal_conductivity]
        type = PorousFlowThermalConductivityIdeal
        dry_thermal_conductivity = '4 0 0  0 4 0  0 0 4' # W/m/K
      []  
  []
  
  [Preconditioning]
    [smp]
      type = SMP
      full = true
      petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
      petsc_options_value = ' lu       mumps'
    []
  []
  
  [Executioner]
    type = Transient
    solve_type = NEWTON
    end_time = 3.0e9
    line_search = none
    dtmin = 100
    [TimeStepper]
      type = IterationAdaptiveDT
      dt = 1000
    []
    [Adaptivity]
        interval = 1
        refine_fraction = 0.2
        coarsen_fraction = 0.3
        max_h_level = 4
      []
  []
  
  [Postprocessors]
   [max_gassat] #check boundary behaves as expected
        type = ElementExtremeValue
        variable = gas_sat
        execute_on = 'initial timestep_end'
    []
[]

[Controls]
    [period0]
        type = TimePeriod
        disable_objects = 'BoundaryCondition::t_dike_neumann'
        enable_objects = 'BoundaryCondition::t_dike_dirichlet'
        start_time = '0'
        end_time = '63072000'
        execute_on = 'initial timestep_begin'
      []
    
      [period2]
        type = TimePeriod
        disable_objects = 'BoundaryCondition::t_dike_dirichlet'
        enable_objects = 'BoundaryCondition::t_dike_neumann'
        start_time = '63072000'
        end_time = '3e9'
        execute_on = 'initial timestep_begin'
      []
[]

  [Outputs]
    perf_graph = true
    exodus = true
    execute_on = 'initial timestep_end failed'
    [residuals]
        type = TopResidualDebugOutput
        num_residuals = 1
        execute_on = 'NONLINEAR TIMESTEP_END INITIAL'
    []
  []
  
  [Debug]
    show_var_residual_norms = true
  []

Thank you for your attention!

[cpgr]

Just letting you know that I've seen this but am away until next week so can't look at it until then. Sounds interesting though so make sure that you bug me if I forget to look at it early next week 😅

[GiudGiud]

So the fix would be for TabulatedFluidProperties to check for those implementations in _fp and use the ones implemented there

this could be fine as a backup method.
a boolean parameter like "use_parameter_fp_for_missing_implementations" and then explicitly forwarding to _fp for all the ones we need when it is set to true

[cpgr]

This is probably my fault! I'm happy to help make this work.

Also, I totally forgot about your earlier question. I'm really sorry about that!

[aikubo]

I opened a new issue about TabulatedFluidProperties: #27457

[aikubo]

This is probably my fault! I'm happy to help make this work.

Also, I totally forgot about your earlier question. I'm really sorry about that!

No worries @cpgr ! I hope that using tabulated fluid properties will help but I'm not positive it will. The blow up of the residual for enthalpy is my major concern which I noted above. Do you have insight to what causes that?

[aikubo]

Adding tabulated fluid properties doesn't seem to help with convergence. Any help on this @cpgr ?

this pull request is mostly done #27506