From 2b9ad77a399718d0b869b4241a454e70e2f7f746 Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Fri, 3 Nov 2023 12:55:35 +0100 Subject: [PATCH] KINSOL: option to select QR orthogonalization strategy This option is available since SUNDIALS 6.0.0. --- include/deal.II/sundials/kinsol.h | 58 ++++++++++++++++++++++++++++- source/sundials/kinsol.cc | 61 ++++++++++++++++++++++++++----- 2 files changed, 108 insertions(+), 11 deletions(-) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 1c980951da..f7f57ceac6 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/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; }; /** diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index d4ec60ab9f..7eb334c159 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -55,16 +55,17 @@ namespace SUNDIALS { template KINSOL::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, -- 2.39.5