Convergence issues after adapting NS AD kernels for incompressible Stokes (PSPG/SUPG enabled)

[bo-qian]

Check these boxes if you have followed the posting rules.

• Q&A Navier Stokes is the most appropriate category for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (avoid screenshots if possible, triple back quotes before/after pasted text, etc)

Question

Note:

This question is based on the previous one (#31590). In the previous case, all kernels related to Stokes and incompressibility were written by me (hand-written Jacobian, non-AD kernels), but I could not resolve some near-zero pressure issues. Therefore, yesterday I copied the relevant official AD kernels and modified parts to implement my model. The main files I copied and modified are the Materials: INSADMaterial.C/.h, INSADTauMaterial.C/.h, the kernel INSADMomentumPressure.C/.h, and I also added a new kernel INSADMomentumSurfaceTension.C/.h. Since my incompressibility equation needs to be multiplied by an interpolation function and the viscosity is variable (depends on C), I changed the strong-form expressions in the Material classes to match what I need.

Used input files:

mu_input = "viscosity_sintering_mu_2D.i"

[Mesh]
[]

[AuxVariables]
  [vel_x]
    order = SECOND
    family = LAGRANGE
  []
  [vel_y]
    order = SECOND
    family = LAGRANGE
  []
  [c]
    order = FIRST
    family = LAGRANGE
  []
  [mu_c]
    order = FIRST
    family = LAGRANGE
  []
  [c_mid]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxKernels]
  [vel_x]
    type = VectorVariableComponentAux
    variable = vel_x
    vector_variable = velocity
    component = 'x'
  []
  [vel_y]
    type = VectorVariableComponentAux
    variable = vel_y
    vector_variable = velocity
    component = 'y'
  []
[]

[Variables]
  [./velocity]
    order = SECOND
    family = LAGRANGE_VEC
  [../]
  [./p]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Kernels]
  [./mass]
    type = INSADMass
    variable = p
  [../]
  # [mass_kernel]
  #   type = MassMatrix
  #   variable = p
  #   matrix_tags = 'mass'
  # []
  [./mass_pspg]
    type = INSADMassPSPG
    variable = p
  [../]

  [./momentum_viscous]
    type = INSADMomentumViscous
    variable = velocity
    viscous_form = traction
  [../]

  [./momentum_pressure]
    type = INSADMomentumPressureDG
    variable = velocity
    pressure = p
    integrate_p_by_parts = true
  [../]

  [./momentum_surface_tension]
    type = INSADMomentumSurfaceTensionDG
    variable = velocity
  [../]

  [./momentum_supg]
    type = INSADMomentumSUPG
    variable = velocity
    velocity = velocity
  [../]
[]

# [Problem]
#   type = NavierStokesProblem
#   mass_matrix = 'mass'
#   extra_tag_matrices = 'mass'
#   use_pressure_mass_matrix = true
# []

[BCs]
  [bcs_velocity]
    type = ADVectorFunctionDirichletBC
    variable = velocity
    boundary = '0 1 2 3'
  []
  [./pressure_pin_1]
    type = DirichletBC
    variable = p
    boundary = 'pinned_node_1'
    value = 0
  [../]
[]

[Materials]
  [./const]
    type = ADGenericConstantMaterial
    prop_names = 'rho'
    prop_values = '1'
  [../]
  [ins_mat]
    type = INSADTauMaterialDG
    velocity = velocity
    pressure = p
    phase_field_mid = c_mid
    chemical_potential = mu_c
  []
[]

[MultiApps]
  [./Mu_Initial]
    type = FullSolveMultiApp
    input_files = ${mu_input}
    execute_on = 'INITIAL'
    clone_parent_mesh = true
  [../]
[]

[Transfers]
  [./MuToStokes]
    type = MultiAppCopyTransfer
    from_multi_app = Mu_Initial
    source_variable = 'c_mid_initial c_mid_initial mu' 
    variable = 'c c_mid mu_c'
    # execute_on = 'INITIAL'
  [../]
[]

[Preconditioning]
  [./pc]
    type = SMP
    full = true
    petsc_options = '-ksp_converged_reason'
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -ksp_gmres_restart -pc_factor_shift_type'
    petsc_options_value = 'lu       mumps                      2500               NONZERO'
  [../]
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
  verbose = true
  
  # automatic_scaling = true

  # nl_rel_tol = 1e-15
  nl_abs_tol = 1e-7
[]

[Outputs]
  exodus = true
  print_linear_residuals = true
[]

[Debug]
  show_var_residual_norms = true
[]

Different solver parameters and the corresponding forms of non-convergence:

When using the lu preconditioning,

[Preconditioning]
  [./pc]
    type = SMP
    full = true
    petsc_options = '-ksp_converged_reason'
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -ksp_gmres_restart -pc_factor_shift_type'
    petsc_options_value = 'lu       mumps                      2500               NONZERO'
  [../]
[]

The output log shows,

Executing MultiApps on INITIAL
Stokes0:     |residual|_2 of individual variables:
Stokes0:                    velocity: 368.69
Stokes0:                    p:        0
Stokes0:  0 Nonlinear |R| = 3.686897e+02
Stokes0:       0 Linear |R| = 3.686897e+02
Stokes0:       1 Linear |R| = 6.255615e-09
Stokes0:     |residual|_2 of individual variables:
Stokes0:                    velocity: 366.004
Stokes0:                    p:        38.3905
Stokes0:  1 Nonlinear |R| = 3.680122e+02
Stokes0:       0 Linear |R| = 3.680122e+02
Stokes0:       1 Linear |R| = 1.256362e-05
Stokes0:     |residual|_2 of individual variables:
Stokes0:                    velocity: 366.004
Stokes0:                    p:        38.3893
Stokes0:  2 Nonlinear |R| = 3.680121e+02
Stokes0:       0 Linear |R| = 3.680121e+02
Stokes0:       1 Linear |R| = 2.582907e-05
Stokes0:     |residual|_2 of individual variables:
Stokes0:                    velocity: 366.025
Stokes0:                    p:        37.822
Stokes0:  3 Nonlinear |R| = 3.679744e+02
Stokes0:       0 Linear |R| = 3.679744e+02
Stokes0:       1 Linear |R| = 5.060723e-07

When using FSP preconditioning and enable MassMatrix , [Problem]/ type = NavierStokesProblem/ mass_matrix = 'mass'.

[Preconditioning]
  [FSP]
    type = FSP
    topsplit = 'up'
    [up]
      splitting = 'vel p'
      splitting_type  = schur
      petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side -ksp_rtol'
      petsc_options_value = 'full                            self                              300                fgmres    right        1e-4'
    []
      [vel]
        vars = 'velocity'
        petsc_options = '-ksp_converged_reason'
        petsc_options_iname = '-pc_type -pc_hypre_type -ksp_type -ksp_rtol -ksp_gmres_restart -ksp_pc_side'
        petsc_options_value = 'hypre    boomeramg      gmres     1e-2      300                right'
      []
      [p]
        vars = 'p'
        petsc_options = '-ksp_converged_reason'
        petsc_options_iname = '-ksp_type -ksp_gmres_restart -ksp_rtol -pc_type -ksp_pc_side -pc_jacobi_type'
        petsc_options_value = 'fgmres    300                1e-2      jacobi   right        rowsum'
      []
  []
[]

The output log shows that it keeps iterating, but there is no sign of the residual decreasing.

Executing MultiApps on INITIAL
Stokes0:     |residual|_2 of individual variables:
Stokes0:                    velocity: 368.69
Stokes0:                    p:        0
Stokes0:  0 Nonlinear |R| = 3.686897e+02
Stokes0:       0 Linear |R| = 3.686897e+02
      Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 5
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 42
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 7
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        Linear fieldsplit_vel_ solve converged due to CONVERGED_RTOL iterations 4
        ······

Based on the above, in what areas can I make improvements, and what methods can I use to enhance convergence to achieve a successful computation?

Additional information

Mesh size and type: 225*225, Taylor-Hood
Reynolds number: None
Discretization (finite element CG/DG, finite volume, etc): finite element CG
Models (turbulence, porous media, etc): Stokes Flow
Solver method (fully coupled, segregated, multiapps, etc): multiapps
Base input you started from: Already given

[GiudGiud]

Hello

ince my incompressibility equation needs to be multiplied by an interpolation function

Did the simulation solve properly when you don't do that?
You need to start from a working input and modify it to include the models you want.

[bo-qian]

Thank you for your response. Indeed, my overall governing equations were somewhat singular. I have now resolved the issue where only single-core computation was possible. The solution was to add -ksp_type preonly to the solver settings at the beginning, which allowed for multi-core parallel computation. This situation is based on my disabling of PSPG and SUPG, as well as using Taylor-Hood elements for the calculation.

[Preconditioning]
  [Direct]
    type = SMP
    full = true
    petsc_options = '-ksp_converged_reason'
    petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_type -pc_factor_shift_type -mat_mumps_icntl_7 -mat_mumps_icntl_24'
    petsc_options_value = 'preonly   lu       mumps                      NONZERO               7                  1'
  []
[]

Can I ask why it works when I add the preonly setting?

[lindsayad]

Preonly will disable a GMRES iterative linear Krylov solver. In general using GMRES should not impact your ability to run in parallel. Do you get an error in parallel when you use GMRES?

[bo-qian]

Previously I was using this setup. Regardless of whether I used a hand-coded Jacobian or automatic differentiation (AD), it would only run on a single core. Running on multiple cores triggered memory errors and the near-zero value error mentioned earlier.

[Preconditioning]
  [CH_Stokes]
    type = SMP
    full = true
    petsc_options = '-ksp_converged_reason'
    petsc_options_iname = '-pc_type -pc_factor_shift_type'
    petsc_options_value = 'lu       NONZERO'
  []
[]

[lindsayad]

So you're saying that if you use

[Preconditioning]
  [Direct]
    type = SMP
    full = true
    petsc_options = '-ksp_converged_reason'
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_type -pc_factor_shift_type -mat_mumps_icntl_7 -mat_mumps_icntl_24'
    petsc_options_value = 'lu       mumps                      NONZERO               7                  1'
  []
[]

then you get something like FACTOR_OUTMEMORY when you run in parallel? If that is the case you might additionally add the option -mat_mumps_icntl_14 300 which would bump your memory working space by 300%. You can look at this PETSc man page for all the mumps options

[bo-qian]

OK, thank you, I will learn more in this page.