Enables libmesh preconditioner reuse capability

        closes idaholab#21868

framework/doc/content/source/systems/NonlinearSystem.md

@@ -323,6 +323,73 @@ only works for symmetric positive-definite matrices. Because of its greater
 flexibility, GMRES is the default linear solution algorithm in PETSc and
 consequently for MOOSE.
 
+## Reusing preconditioners id=reuse_preconditioners
+
+The simple version of GRMES and other iterative methods converge only very
+slowly.  To improve convergence PETSc and other iterative solver packages
+apply a [preconditioner](https://en.wikipedia.org/wiki/Preconditioner) to the
+system of equations/sparse matrix before applying the iterative solver.
+
+A great number of preconditioners exist, but 
+[multigrid](https://en.wikipedia.org/wiki/Multigrid_method)
+methods are often among the best choices.  
+The [HYPRE](application_development/hypre.md optional=true) package, 
+specifically the 
+BoomerAMG preconditioner, is often a good choice for a preconditioner to
+condition the system of equations resulting from the MOOSE simulation.
+
+A direct factorization of the sparse system of equations is a *very*
+good preconditioner.  An iterative method using the factorized matrix
+as a preconditioner will typically converge to machine precision in a single
+iteration.  However, as noted above, factorizing the sparse system of equations
+for a large simulation is numerically expensive.
+
+One option is to reuse a preconditioner for many iterations of the 
+linear iterative solver.  The preconditioner can be carried over
+through nonlinear iterations and even across load steps.  MOOSE
+allows the user to do this with the `reuse_preconditioner` flag.
+Setting
+
+```
+  reuse_preconditioner = true
+  reuse_preconditioner_max_its = 20
+```
+
+in the `[Executioner]` block will reuse the same preconditioner until
+the number of linear iterations required to solve the linearized system of
+equations exceeds 20.  At that time the solution algorithm will reform
+the preconditioner and repeat the process.
+
+Using these parameters in combination with a direct factorization of the
+system can be very efficient.  The following is an example of how to
+configure PETSc and MOOSE to solve the equations with this combination:
+
+```
+  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_type'
+  petsc_options_value = 'lu superlu_dist gmres'
+  l_tol = 1e-8
+  l_max_its = 100
+
+  reuse_preconditioner = false
+  reuse_preconditioner_max_its = 20
+
+  nl_max_its = 10
+  nl_rel_tol = 1e-8
+  nl_abs_tol = 1e-10
+```
+
+This solver strategy can be very efficient when solving a problem on a 
+desktop or workstation computer (i.e. not on a cluster).  However, it
+is a heuristic -- the efficiency of the approach and the optimal value
+of `reuse_preconditioner_max_its` will vary greatly depending on the type of
+physics you are solving and on how quickly the structure of the problem
+changes from load step to load step.
+
+On larger machines reusing the AMG preconditioner provided by BoomerAMG
+may be a way to increase the efficiency of the overall simulation.  Again,
+this is a heuristic and the effectiveness of the approach and good values
+of `reuse_preconditioner_max_its` are going to vary.
+
 ## Augmenting Sparsity id=augmenting_sparsity
 
 One such routine is `NonlinearSystemBase::augmentSparsity`, which as its name
@@ -339,29 +406,3 @@ information comes from [`NearestNodeLocators`](/NearestNodeLocator.md) held by
 displacements). In the future sparsity augmentation from constraints will occur
 through [`RelationshipManagers`](/RelationshipManager.md) instead of through the
 `augmentSparsity` method.
-
-## Computing Residual and Jacobian Together id=resid_and_jac_together
-
-The default behavior in MOOSE is to have separate functions compute the residual
-and Jacobian. However, with the advent of [NonlinearSystem.md#AD] it can make
-sense to use a single function to compute the residual and Jacobian
-simultaneously. At the local residual object level, automatic differentiation (AD)
-already computes the residual and Jacobian simultaneously, with the dual number
-at the core of AD holding value (residual) and derivatives
-(Jacobian). Simultaneous evaluation of residual and Jacobian using a single
-function can be triggered by setting
-[!param](/Executioner/Steady/residual_and_jacobian_together) to `true`. What this does in
-the background is funnels the (generally AD) computed local residuals and
-Jacobians into the global residual vector and Jacobian
-matrix respectively when PETSc calls the libMesh/MOOSE residual/function
-evaluation routine. Then when PETSc calls the libMesh/MOOSE Jacobian evaluation
-routine, we simply return because the global matrix has already been computed.
-
-Computing the residual and Jacobian together has shown 20% gains for
-Navier-Stokes finite volume simulations for which automatic differentiation is
-leveraged even during standard residual evaluations. Other areas where computing
-the residual and Jacobian together may be advantageous is in simulations in
-which there are quite a few nonlinear iterations per timestep, for which the
-cost of an additional Jacobian evaluation during the final residual evaluation
-is amortized. Also, simulations in which material property calculations are very
-expensive may be good candidates for computing the residual and Jacobian together.

framework/src/executioners/FEProblemSolve.C

@@ -126,11 +126,25 @@ FEProblemSolve::validParams()
                         false,
                         "Whether to compute the residual and Jacobian together.");
 
+  params.addParam<bool>("reuse_preconditioner",
+                        false,
+                        "If true reuse the previously calculated "
+                        "preconditioner for the linearized "
+                        "system across multiple solvers, "
+                        "controlled by reuse_preconditioner_max_its");
+  params.addParam<unsigned int>("reuse_preconditioner_max_its",
+                                25,
+                                "Reuse the previously calculated "
+                                "preconditioner for the linear system "
+                                "until the number of linear iterations "
+                                "exceeds this number");
+
   params.addParamNamesToGroup("solve_type l_tol l_abs_tol l_max_its nl_max_its nl_max_funcs "
                               "nl_abs_tol nl_rel_tol nl_abs_step_tol nl_rel_step_tol "
                               "snesmf_reuse_base compute_initial_residual_before_preset_bcs "
                               "num_grids nl_div_tol nl_abs_div_tol residual_and_jacobian_together "
-                              "n_max_nonlinear_pingpong",
+                              "n_max_nonlinear_pingpong reuse_preconditioner "
+                              "reuse_preconditioner_max_its",
                               "Solver");
   params.addParamNamesToGroup(
       "automatic_scaling compute_scaling_once off_diagonals_in_auto_scaling "
@@ -183,6 +197,11 @@ FEProblemSolve::FEProblemSolve(Executioner & ex)
   es.parameters.set<Real>("nonlinear solver relative step tolerance") =
       getParam<Real>("nl_rel_step_tol");
 
+  es.parameters.set<bool>("reuse preconditioner") = getParam<bool>("reuse_preconditioner");
+
+  es.parameters.set<unsigned int>("reuse preconditioner maximum iterations") =
+      getParam<unsigned int>("reuse_preconditioner_max_its");
+
   _nl._compute_initial_residual_before_preset_bcs =
       getParam<bool>("compute_initial_residual_before_preset_bcs");
 

modules/tensor_mechanics/test/tests/preconditioner_reuse/convergence.i

@@ -0,0 +1,164 @@
+# Simple 3D test
+
+[GlobalParams]
+  displacements = 'disp_x disp_y disp_z'
+  large_kinematics = true
+[]
+
+[Variables]
+  [disp_x]
+  []
+  [disp_y]
+  []
+  [disp_z]
+  []
+[]
+
+[Mesh]
+  [msh]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+[]
+
+[Kernels]
+  [sdx]
+    type = TotalLagrangianStressDivergence
+    variable = disp_x
+    component = 0
+  []
+  [sdy]
+    type = TotalLagrangianStressDivergence
+    variable = disp_y
+    component = 1
+  []
+  [sdz]
+    type = TotalLagrangianStressDivergence
+    variable = disp_z
+    component = 2
+  []
+[]
+
+[Functions]
+  [pullx]
+    type = ParsedFunction
+    value = '4000 * t'
+  []
+  [pully]
+    type = ParsedFunction
+    value = '-2000 * t'
+  []
+  [pullz]
+    type = ParsedFunction
+    value = '3000 * t'
+  []
+  [lambda_function]
+    type = ParsedFunction
+    value = '1000.0*(t+1.0)'
+  []
+[]
+
+[BCs]
+  [leftx]
+    type = DirichletBC
+    preset = true
+    boundary = left
+    variable = disp_x
+    value = 0.0
+  []
+  [lefty]
+    type = DirichletBC
+    preset = true
+    boundary = left
+    variable = disp_y
+    value = 0.0
+  []
+  [leftz]
+    type = DirichletBC
+    preset = true
+    boundary = left
+    variable = disp_z
+    value = 0.0
+  []
+  [pull_x]
+    type = FunctionNeumannBC
+    boundary = right
+    variable = disp_x
+    function = pullx
+  []
+  [pull_y]
+    type = FunctionNeumannBC
+    boundary = top
+    variable = disp_y
+    function = pully
+  []
+  [pull_z]
+    type = FunctionNeumannBC
+    boundary = right
+    variable = disp_z
+    function = pullz
+  []
+[]
+
+[Materials]
+  [compute_stress]
+    type = ComputeNeoHookeanStress
+    lambda = lambda
+    mu = 67000.0
+  []
+  [lambda]
+    type = GenericFunctionMaterial
+    prop_names = lambda
+    prop_values = lambda_function
+  []
+  [compute_strain]
+    type = ComputeLagrangianStrain
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+
+  solve_type = 'newton'
+  line_search = none
+
+  petsc_options = ''
+  petsc_options_iname = '-pc_type -ksp_type'
+  petsc_options_value = 'lu gmres'
+  l_tol = 1e-8
+  l_max_its = 100
+
+  reuse_preconditioner = false
+  reuse_preconditioner_max_its = 20
+
+  nl_max_its = 10
+  nl_rel_tol = 1e-8
+  nl_abs_tol = 1e-10
+
+  start_time = 0.0
+  dt = 1.0
+  dtmin = 1.0
+  end_time = 10.0
+[]
+
+[Reporters/iteration_info]
+  type = IterationInfo
+[]
+
+[Outputs]
+  exodus = false
+  [./csv]
+    type = CSV
+    file_base = base_case
+  [../]
+[]

modules/tensor_mechanics/test/tests/preconditioner_reuse/gold/base_case.csv

@@ -0,0 +1,12 @@
+time,iteration_info/num_fixed_point_iterations,iteration_info/num_linear_iterations,iteration_info/num_nonlinear_iterations,iteration_info/time,iteration_info/timestep
+0,0,0,0,0,0
+1,1,4,4,1,1
+2,1,4,4,2,2
+3,1,4,4,3,3
+4,1,4,4,4,4
+5,1,3,3,5,5
+6,1,3,3,6,6
+7,1,3,3,7,7
+8,1,3,3,8,8
+9,1,3,3,9,9
+10,1,3,3,10,10

modules/tensor_mechanics/test/tests/preconditioner_reuse/gold/reuse_case.csv

@@ -0,0 +1,12 @@
+time,iteration_info/num_fixed_point_iterations,iteration_info/num_linear_iterations,iteration_info/num_nonlinear_iterations,iteration_info/time,iteration_info/timestep
+0,0,0,0,0,0
+1,1,29,4,1,1
+2,1,47,4,2,2
+3,1,56,4,3,3
+4,1,64,4,4,4
+5,1,52,3,5,5
+6,1,55,3,6,6
+7,1,58,3,7,7
+8,1,42,3,8,8
+9,1,16,3,9,9
+10,1,20,3,10,10

modules/tensor_mechanics/test/tests/preconditioner_reuse/tests

@@ -0,0 +1,19 @@
+[Tests]
+  issues = '#21868'
+  [without_reuse]
+    type = CSVDiff
+    input = convergence.i
+    csvdiff = base_case.csv
+    requirement = 'Convergence matches previous version of MOOSE without the '
+                  'preconditioner reuse system'
+  []
+  [with_reuse]
+    type = CSVDiff
+    input = convergence.i
+    cli_args = "Executioner/reuse_preconditioner=true Outputs/csv/file_base=reuse_case"
+    csvdiff = reuse_case.csv
+    requirement = 'Preconditioner is reused until the linear iterations exceed '
+                  'the value of reuse_preconditioner_max_its upon which the '
+                  'system recalculates the preconditioner'
+  []
+[]