Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines DataBox tags for the linear solver
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <cstddef>
      10             : #include <string>
      11             : 
      12             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      13             : #include "DataStructures/DataBox/Tag.hpp"
      14             : #include "Options/String.hpp"
      15             : #include "Utilities/Gsl.hpp"
      16             : #include "Utilities/PrettyType.hpp"
      17             : 
      18             : /// Functionality for solving nonlinear systems of equations
      19             : namespace NonlinearSolver {
      20             : 
      21             : /// Options related to nonlinear solvers
      22           1 : namespace OptionTags {
      23             : 
      24             : /*!
      25             :  * \brief Sufficient decrease parameter of the line search globalization
      26             :  *
      27             :  * The sufficient decrease parameter is the acceptable decrease of the residual
      28             :  * magnitude in each step of the nonlinear solver. It is measured as a fraction
      29             :  * of the predicted decrease in residual magnitude if the problem was linear.
      30             :  * For example, a sufficient decrease parameter of 1 means that a nonlinear
      31             :  * solver step is expected to decrease the residual exactly as expected for a
      32             :  * linear problem, i.e. immediately to zero. A sufficient decrease parameter of
      33             :  * 0.5 means that decreasing the residual by half of that amount in each
      34             :  * nonlinear solver step is acceptable.
      35             :  *
      36             :  * Nonlinear solver steps that fail the sufficient decrease condition (also
      37             :  * known as _Armijo condition_) undergo a globalization procedure such as a line
      38             :  * search.
      39             :  *
      40             :  * A typical value for the sufficient decrease parameter is \f$10^{-4}\f$. Set
      41             :  * to values closer to unity when the nonlinear solver overshoots, e.g. when the
      42             :  * initial guess is particularly bad. Larger values mean the nonlinear solver is
      43             :  * stricter with accepting steps, preferring to apply the globalization
      44             :  * strategy.
      45             :  */
      46             : template <typename OptionsGroup>
      47           1 : struct SufficientDecrease {
      48           0 :   using type = double;
      49           0 :   static constexpr Options::String help = {
      50             :       "Fraction of decrease predicted by linearization"};
      51           0 :   static type lower_bound() { return 0.; }
      52           0 :   static type upper_bound() { return 1.; }
      53           0 :   static type suggested_value() { return 1.e-4; }
      54           0 :   using group = OptionsGroup;
      55             : };
      56             : 
      57             : /*!
      58             :  * \brief Nonlinear solver steps are damped by this factor
      59             :  *
      60             :  * Instead of attempting to take full-length steps when correcting the solution
      61             :  * in each nonlinear solver step (see `NonlinearSolver::Tags::Correction`),
      62             :  * reduce the step length by this factor. This damping occurs before any
      63             :  * globalization steps that may further reduce the step length.
      64             :  */
      65             : template <typename OptionsGroup>
      66           1 : struct DampingFactor {
      67           0 :   using type = double;
      68           0 :   static constexpr Options::String help = {
      69             :       "Multiply corrections by this factor"};
      70           0 :   static type lower_bound() { return 0.; }
      71           0 :   static type upper_bound() { return 1.; }
      72           0 :   static type suggested_value() { return 1.; }
      73           0 :   using group = OptionsGroup;
      74             : };
      75             : 
      76             : /*!
      77             :  * \brief The maximum number of allowed globalization steps
      78             :  *
      79             :  * Nonlinear solves of well-posed problems should never hit this limit because
      80             :  * the step size shrinks and eventually triggers the sufficient-decrease
      81             :  * condition (see `NonlinearSolver::OptionTags::SufficientDecrease`). So the
      82             :  * suggested value just provides a safety-net to prevent the globalization
      83             :  * from running forever when the problem is ill-posed.
      84             :  */
      85             : template <typename OptionsGroup>
      86           1 : struct MaxGlobalizationSteps {
      87           0 :   using type = size_t;
      88           0 :   static constexpr Options::String help = {
      89             :       "Maximum number of globalization steps"};
      90           0 :   static type suggested_value() { return 40; }
      91           0 :   using group = OptionsGroup;
      92             : };
      93             : 
      94             : }  // namespace OptionTags
      95             : 
      96           0 : namespace Tags {
      97             : 
      98             : /*!
      99             :  * \brief The correction \f$\delta x\f$ to improve a solution \f$x_0\f$
     100             :  *
     101             :  * A linear problem \f$Ax=b\f$ can be equivalently formulated as the problem
     102             :  * \f$A\delta x=b-A x_0\f$ for the correction \f$\delta x\f$ to an initial guess
     103             :  * \f$x_0\f$. More importantly, we can use a correction scheme to solve a
     104             :  * nonlinear problem \f$A_\mathrm{nonlinear}(x)=b\f$ by repeatedly solving a
     105             :  * linearization of it. For instance, a Newton-Raphson scheme iteratively
     106             :  * refines an initial guess \f$x_0\f$ by repeatedly solving the linearized
     107             :  * problem
     108             :  *
     109             :  * \f{equation}
     110             :  * \frac{\delta A_\mathrm{nonlinear}}{\delta x}(x_k)\delta x_k =
     111             :  * b-A_\mathrm{nonlinear}(x_k)
     112             :  * \f}
     113             :  *
     114             :  * for the correction \f$\delta x_k\f$ and then updating the solution as
     115             :  * \f$x_{k+1}=x_k + \delta x_k\f$.
     116             :  */
     117             : template <typename Tag>
     118           1 : struct Correction : db::PrefixTag, db::SimpleTag {
     119           0 :   using type = typename Tag::type;
     120           0 :   using tag = Tag;
     121             : };
     122             : 
     123             : /*!
     124             :  * \brief The nonlinear operator \f$A_\mathrm{nonlinear}\f$ applied to the data
     125             :  * in `Tag`
     126             :  */
     127             : template <typename Tag>
     128           1 : struct OperatorAppliedTo : db::PrefixTag, db::SimpleTag {
     129           0 :   static std::string name() {
     130             :     // Add "Nonlinear" prefix to abbreviate the namespace for uniqueness
     131             :     return "NonlinearOperatorAppliedTo(" + db::tag_name<Tag>() + ")";
     132             :   }
     133           0 :   using type = typename Tag::type;
     134           0 :   using tag = Tag;
     135             : };
     136             : 
     137             : /*!
     138             :  * \brief The nonlinear residual
     139             :  * \f$r_\mathrm{nonlinear} = b - A_\mathrm{nonlinear}(\delta x)\f$
     140             :  */
     141             : template <typename Tag>
     142           1 : struct Residual : db::PrefixTag, db::SimpleTag {
     143           0 :   static std::string name() {
     144             :     // Add "Nonlinear" prefix to abbreviate the namespace for uniqueness
     145             :     return "NonlinearResidual(" + db::tag_name<Tag>() + ")";
     146             :   }
     147           0 :   using type = typename Tag::type;
     148           0 :   using tag = Tag;
     149             : };
     150             : 
     151             : /// Compute the residual \f$r=b - Ax\f$ from the `SourceTag` \f$b\f$ and the
     152             : /// `db::add_tag_prefix<NonlinearSolver::Tags::OperatorAppliedTo, FieldsTag>`
     153             : /// \f$Ax\f$.
     154             : template <typename FieldsTag, typename SourceTag>
     155           1 : struct ResidualCompute : db::add_tag_prefix<Residual, FieldsTag>,
     156             :                          db::ComputeTag {
     157           0 :   using base = db::add_tag_prefix<Residual, FieldsTag>;
     158           0 :   using argument_tags =
     159             :       tmpl::list<SourceTag, db::add_tag_prefix<OperatorAppliedTo, FieldsTag>>;
     160           0 :   using return_type = typename base::type;
     161           0 :   static void function(
     162             :       const gsl::not_null<return_type*> residual,
     163             :       const typename SourceTag::type& source,
     164             :       const typename db::add_tag_prefix<OperatorAppliedTo, FieldsTag>::type&
     165             :           operator_applied_to_fields) {
     166             :     *residual = source - operator_applied_to_fields;
     167             :   }
     168             : };
     169             : 
     170             : /*!
     171             :  * \brief The length of nonlinear solver steps
     172             :  *
     173             :  * Instead of taking full-length nonlinear solver steps when correcting the
     174             :  * solution as detailed in `NonlinearSolver::Tags::Correction`, the correction
     175             :  * is multiplied by this step length.
     176             :  *
     177             :  * The `NonlinearSolver::Tags::DampingFactor` multiplies the initial length of
     178             :  * each step such that the nonlinear solver never takes full-length steps if the
     179             :  * damping factor is below one. The step length can be further reduced by the
     180             :  * globalization procedure. See `NonlinearSolver::NewtonRaphson` for details.
     181             :  */
     182             : template <typename OptionsGroup>
     183           1 : struct StepLength : db::SimpleTag {
     184           0 :   static std::string name() {
     185             :     return "StepLength(" + pretty_type::name<OptionsGroup>() + ")";
     186             :   }
     187           0 :   using type = double;
     188             : };
     189             : 
     190             : /*!
     191             :  * \brief Sufficient decrease parameter of the line search globalization
     192             :  *
     193             :  * \see `NonlinearSolver::OptionTags::SufficientDecrease`
     194             :  */
     195             : template <typename OptionsGroup>
     196           1 : struct SufficientDecrease : db::SimpleTag {
     197           0 :   static std::string name() {
     198             :     return "SufficientDecrease(" + pretty_type::name<OptionsGroup>() + ")";
     199             :   }
     200           0 :   using type = double;
     201           0 :   static constexpr bool pass_metavariables = false;
     202           0 :   using option_tags = tmpl::list<OptionTags::SufficientDecrease<OptionsGroup>>;
     203           0 :   static type create_from_options(const type& option) { return option; }
     204             : };
     205             : 
     206             : /*!
     207             :  * \brief Nonlinear solver steps are damped by this factor
     208             :  *
     209             :  * \see `NonlinearSolver::OptionTags::DampingFactor`
     210             :  */
     211             : template <typename OptionsGroup>
     212           1 : struct DampingFactor : db::SimpleTag {
     213           0 :   static std::string name() {
     214             :     return "DampingFactor(" + pretty_type::name<OptionsGroup>() + ")";
     215             :   }
     216           0 :   using type = double;
     217           0 :   static constexpr bool pass_metavariables = false;
     218           0 :   using option_tags = tmpl::list<OptionTags::DampingFactor<OptionsGroup>>;
     219           0 :   static type create_from_options(const type& option) { return option; }
     220             : };
     221             : 
     222             : /*!
     223             :  * \brief The maximum number of allowed globalization steps
     224             :  *
     225             :  * \see `NonlinearSolver::OptionTags::MinStepLength`
     226             :  */
     227             : template <typename OptionsGroup>
     228           1 : struct MaxGlobalizationSteps : db::SimpleTag {
     229           0 :   static std::string name() {
     230             :     return "MaxGlobalizationSteps(" + pretty_type::name<OptionsGroup>() + ")";
     231             :   }
     232           0 :   using type = size_t;
     233           0 :   static constexpr bool pass_metavariables = false;
     234           0 :   using option_tags =
     235             :       tmpl::list<OptionTags::MaxGlobalizationSteps<OptionsGroup>>;
     236           0 :   static type create_from_options(const type& option) { return option; }
     237             : };
     238             : 
     239             : /// Prefix indicating the `Tag` is related to the globalization procedure
     240             : template <typename Tag>
     241           1 : struct Globalization : db::PrefixTag, db::SimpleTag {
     242           0 :   using type = typename Tag::type;
     243           0 :   using tag = Tag;
     244             : };
     245             : 
     246             : }  // namespace Tags
     247             : }  // namespace NonlinearSolver