KINSOL: option to select QR orthogonalization strategy

include/deal.II/sundials/kinsol.h

@@ -214,6 +214,52 @@ namespace SUNDIALS
         picard = KIN_PICARD,
       };
 
+
+      /**
+       * Orthogonalization strategy used for fixed-point iteration with
+       * Anderson acceleration. While `modified_gram_schmidt` is a good default
+       * choice, it suffers from excessive communication for a growing number of
+       * acceleration vectors. Instead, the user may use one of the other
+       * strategies.  For a detailed documentation consult the SUNDIALS manual.
+       *
+       * @note Similar strategies are defined in
+       * LinearAlgebra::OrthogonalizationStrategy.
+       *
+       * @note The values assigned to the enum constants matter since they are
+       * directly passed to KINSOL.
+       */
+      enum OrthogonalizationStrategy
+      {
+        /**
+         * The default, stable strategy. Requires an increasing number of
+         * communication steps with increasing Anderson subspace size.
+         */
+        modified_gram_schmidt = 0,
+
+        /**
+         * Use a compact WY representation in the orthogonalization process of
+         * the inverse iteration with the Householder transformation. This
+         * strategy requires two communication steps and a very small linear
+         * system solve.
+         */
+        inverse_compact = 1,
+
+        /**
+         * Classical Gram Schmidt, that can work on multi-vectors in parallel
+         * and thus always requires three communication steps. However, it
+         * is less stable in terms of roundoff error propagation, requiring
+         * additional re-orthogonalization steps more frequently.
+         */
+        classical_gram_schmidt = 2,
+
+        /**
+         * Classical Gram Schmidt with delayed re-orthogonalization, which
+         * reduces the number of communication steps from three to two compared
+         * to `classical_gram_schmidt`.
+         */
+        delayed_classical_gram_schmidt = 3,
+      };
+
       /**
        * Initialization parameters for KINSOL.
        *
@@ -241,6 +287,8 @@ namespace SUNDIALS
        * Fixed point and Picard parameters:
        *
        * @param anderson_subspace_size Anderson acceleration subspace size
+       * @param anderson_qr_orthogonalization Orthogonalization strategy for QR
+       * factorization when using Anderson acceleration
        */
       AdditionalData(const SolutionStrategy &strategy = linesearch,
                      const unsigned int maximum_non_linear_iterations = 200,
@@ -251,7 +299,9 @@ namespace SUNDIALS
                      const double       maximum_newton_step           = 0.0,
                      const double       dq_relative_error             = 0.0,
                      const unsigned int maximum_beta_failures         = 0,
-                     const unsigned int anderson_subspace_size        = 0);
+                     const unsigned int anderson_subspace_size        = 0,
+                     const OrthogonalizationStrategy
+                       anderson_qr_orthogonalization = modified_gram_schmidt);
 
       /**
        * Add all AdditionalData() parameters to the given ParameterHandler
@@ -374,6 +424,12 @@ namespace SUNDIALS
        * If you set this to 0, no acceleration is used.
        */
       unsigned int anderson_subspace_size;
+
+      /**
+       * In case Anderson acceleration is used, this parameter selects the
+       * orthogonalization strategy used in the QR factorization.
+       */
+      OrthogonalizationStrategy anderson_qr_orthogonalization;
     };
 
     /**

source/sundials/kinsol.cc

@@ -55,16 +55,17 @@ namespace SUNDIALS
 {
   template <typename VectorType>
   KINSOL<VectorType>::AdditionalData::AdditionalData(
-    const SolutionStrategy &strategy,
-    const unsigned int      maximum_non_linear_iterations,
-    const double            function_tolerance,
-    const double            step_tolerance,
-    const bool              no_init_setup,
-    const unsigned int      maximum_setup_calls,
-    const double            maximum_newton_step,
-    const double            dq_relative_error,
-    const unsigned int      maximum_beta_failures,
-    const unsigned int      anderson_subspace_size)
+    const SolutionStrategy         &strategy,
+    const unsigned int              maximum_non_linear_iterations,
+    const double                    function_tolerance,
+    const double                    step_tolerance,
+    const bool                      no_init_setup,
+    const unsigned int              maximum_setup_calls,
+    const double                    maximum_newton_step,
+    const double                    dq_relative_error,
+    const unsigned int              maximum_beta_failures,
+    const unsigned int              anderson_subspace_size,
+    const OrthogonalizationStrategy anderson_qr_orthogonalization)
     : strategy(strategy)
     , maximum_non_linear_iterations(maximum_non_linear_iterations)
     , function_tolerance(function_tolerance)
@@ -75,6 +76,7 @@ namespace SUNDIALS
     , dq_relative_error(dq_relative_error)
     , maximum_beta_failures(maximum_beta_failures)
     , anderson_subspace_size(anderson_subspace_size)
+    , anderson_qr_orthogonalization(anderson_qr_orthogonalization)
   {}
 
 
@@ -125,6 +127,30 @@ namespace SUNDIALS
     prm.enter_subsection("Fixed point and Picard parameters");
     prm.add_parameter("Anderson acceleration subspace size",
                       anderson_subspace_size);
+
+    static std::string orthogonalization_str("modified_gram_schmidt");
+    prm.add_parameter(
+      "Anderson QR orthogonalization",
+      orthogonalization_str,
+      "Choose among modified_gram_schmidt|inverse_compact|"
+      "classical_gram_schmidt|delayed_classical_gram_schmidt",
+      Patterns::Selection(
+        "modified_gram_schmidt|inverse_compact|classical_gram_schmidt|"
+        "delayed_classical_gram_schmidt"));
+    prm.add_action("Anderson QR orthogonalization",
+                   [&](const std::string &value) {
+                     if (value == "modified_gram_schmidt")
+                       anderson_qr_orthogonalization = modified_gram_schmidt;
+                     else if (value == "inverse_compact")
+                       anderson_qr_orthogonalization = inverse_compact;
+                     else if (value == "classical_gram_schmidt")
+                       anderson_qr_orthogonalization = classical_gram_schmidt;
+                     else if (value == "delayed_classical_gram_schmidt")
+                       anderson_qr_orthogonalization =
+                         delayed_classical_gram_schmidt;
+                     else
+                       Assert(false, ExcInternalError());
+                   });
     prm.leave_subsection();
   }
 
@@ -258,6 +284,21 @@ namespace SUNDIALS
     status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
     AssertKINSOL(status);
 
+#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+    // From the manual: this must be called BEFORE KINInit
+    status = KINSetOrthAA(kinsol_mem, data.anderson_qr_orthogonalization);
+    AssertKINSOL(status);
+#  else
+    AssertThrow(
+      data.anderson_qr_orthogonalization ==
+        AdditionalData::modified_gram_schmidt,
+      ExcMessage(
+        "You specified an orthogonalization strategy for QR factorization "
+        "different from the default (modified Gram-Schmidt) but the installed "
+        "SUNDIALS version does not support this feature. Either choose the "
+        "default or install a SUNDIALS version >= 6.0.0."));
+#  endif
+
     if (data.strategy == AdditionalData::fixed_point)
       status = KINInit(
         kinsol_mem,