Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include "DataStructures/DataBox/PrefixHelpers.hpp"
       7              : #include "DataStructures/DataBox/Prefixes.hpp"
       8              : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/ElementActions.hpp"
       9              : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/ResidualMonitor.hpp"
      10              : #include "ParallelAlgorithms/NonlinearSolver/Tags.hpp"
      11              : #include "Utilities/TMPL.hpp"
      12              : 
      13              : namespace NonlinearSolver {
      14              : /// Items related to the NewtonRaphson nonlinear solver
      15              : namespace newton_raphson {
      16              : 
      17              : /*!
      18              :  * \brief A Newton-Raphson correction scheme for nonlinear systems of equations
      19              :  * \f$A_\mathrm{nonlinear}(x)=b\f$.
      20              :  *
      21              :  * We can use a correction scheme to solve a nonlinear problem
      22              :  * \f$A_\mathrm{nonlinear}(x)=b\f$ by repeatedly solving a linearization of it.
      23              :  * A Newton-Raphson scheme iteratively refines an initial guess \f$x_0\f$ by
      24              :  * repeatedly solving the linearized problem
      25              :  *
      26              :  * \f{equation}
      27              :  * \frac{\delta A_\mathrm{nonlinear}}{\delta x}(x_k) \delta x_k =
      28              :  * b-A_\mathrm{nonlinear}(x_k) \equiv r_k
      29              :  * \f}
      30              :  *
      31              :  * for the correction \f$\delta x_k\f$ and then updating the solution as
      32              :  * \f$x_{k+1}=x_k + \delta x_k\f$.
      33              :  *
      34              :  * The operations we need to supply to the algorithm are the nonlinear operator
      35              :  * \f$A_\mathrm{nonlinear}(x)\f$ and a linear solver for the linearized problem
      36              :  * \f$\frac{\delta A_\mathrm{nonlinear}}{\delta x}(x_k) \delta x_k = r_k\f$.
      37              :  * Each step of the algorithm expects that \f$A_\mathrm{nonlinear}(x)\f$ is
      38              :  * computed and stored in the DataBox as
      39              :  * `db::add_tag_prefix<NonlinearSolver::Tags::OperatorAppliedTo, FieldsTag>`.
      40              :  * To perform a solve, add the `solve` action list to an array parallel
      41              :  * component. Pass the actions that compute \f$A_\mathrm{nonlinear}(x)\f$ as
      42              :  * the first template parameter to `solve`. As the second template parameter,
      43              :  * pass the action list that performs a linear solve of the linearized operator
      44              :  * \f$\frac{\delta A_\mathrm{nonlinear}}{\delta x}(x)\f$ for the field
      45              :  * \f$\delta x\f$ (the `linear_solver_fields_tag`) sourced by
      46              :  * \f$r\f$ (the `linear_solver_source_tag`). You will find suitable iterative
      47              :  * linear solvers in the `LinearSolver` namespace. The third template parameter
      48              :  * allows you to pass any further actions you wish to run in each step of the
      49              :  * algorithm (such as observations). If you add the `solve` action list multiple
      50              :  * times, use the fourth template parameter to label each solve with a different
      51              :  * type.
      52              :  *
      53              :  * \par Globalization:
      54              :  * This nonlinear solver supports a line-search (or "backtracking")
      55              :  * globalization. If a step does not sufficiently decrease the residual (see
      56              :  * `NonlinearSolver::OptionTags::SufficientDecrease` for details on the
      57              :  * sufficient decrease condition), the step length is reduced until the residual
      58              :  * sufficiently decreases or the maximum number of globalization steps is
      59              :  * reached. The reduction in step length is determined by the minimum of a
      60              :  * quadratic (first globalization step) or cubic (subsequent globalization
      61              :  * steps) polynomial interpolation according to Algorithm A6.1.3 in
      62              :  * \cite DennisSchnabel (p. 325) (see
      63              :  * `NonlinearSolver::newton_raphson::next_step_length`). Alternatives to a
      64              :  * line-search globalization, such as a trust-region globalization or more
      65              :  * sophisticated nonlinear preconditioning techniques (see e.g. \cite Brune2015
      66              :  * for an overview), are not currently implemented.
      67              :  *
      68              :  * \par Array sections
      69              :  * This nonlinear solver supports running over a subset of the elements in the
      70              :  * array parallel component (see `Parallel::Section`). Set the
      71              :  * `ArraySectionIdTag` template parameter to restrict the solver to elements in
      72              :  * that section. Only a single section must be associated with the
      73              :  * `ArraySectionIdTag`. The default is `void`, which means running over all
      74              :  * elements in the array. Note that the actions in the `SolveLinearization`
      75              :  * list passed to `solve` will _not_ be restricted to run only on section
      76              :  * elements, so all elements in the array may participate in preconditioning
      77              :  * (see LinearSolver::multigrid::Multigrid).
      78              :  */
      79              : template <typename Metavariables, typename FieldsTag, typename OptionsGroup,
      80              :           typename SourceTag =
      81              :               db::add_tag_prefix<::Tags::FixedSource, FieldsTag>,
      82              :           typename ArraySectionIdTag = void>
      83            1 : struct NewtonRaphson {
      84            0 :   using fields_tag = FieldsTag;
      85            0 :   using options_group = OptionsGroup;
      86            0 :   using source_tag = SourceTag;
      87              : 
      88            0 :   using operand_tag = fields_tag;
      89            0 :   using linear_solver_fields_tag =
      90              :       db::add_tag_prefix<NonlinearSolver::Tags::Correction, fields_tag>;
      91            0 :   using linear_solver_source_tag =
      92              :       db::add_tag_prefix<NonlinearSolver::Tags::Residual, fields_tag>;
      93              : 
      94            0 :   using component_list = tmpl::list<
      95              :       detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
      96              : 
      97            0 :   using initialize_element =
      98              :       detail::InitializeElement<FieldsTag, OptionsGroup, SourceTag>;
      99              : 
     100            0 :   using register_element = tmpl::list<>;
     101              : 
     102            0 :   using amr_projectors = initialize_element;
     103              : 
     104              :   template <typename ApplyNonlinearOperator, typename SolveLinearization,
     105              :             typename ObserveActions = tmpl::list<>,
     106              :             typename Label = OptionsGroup>
     107            0 :   using solve = tmpl::list<
     108              :       detail::PrepareSolve<FieldsTag, OptionsGroup, Label, ArraySectionIdTag>,
     109              :       ApplyNonlinearOperator, ObserveActions,
     110              :       detail::SendInitialResidualMagnitude<FieldsTag, OptionsGroup, Label,
     111              :                                            ArraySectionIdTag>,
     112              :       detail::ReceiveInitialHasConverged<FieldsTag, OptionsGroup, Label,
     113              :                                          ArraySectionIdTag>,
     114              :       detail::PrepareStep<FieldsTag, OptionsGroup, Label, ArraySectionIdTag>,
     115              :       SolveLinearization,
     116              :       detail::PerformStep<FieldsTag, OptionsGroup, Label, ArraySectionIdTag>,
     117              :       ApplyNonlinearOperator,
     118              :       detail::ContributeToResidualMagnitudeReduction<FieldsTag, OptionsGroup,
     119              :                                                      Label, ArraySectionIdTag>,
     120              :       detail::Globalize<FieldsTag, OptionsGroup, Label, ArraySectionIdTag>,
     121              :       ObserveActions,
     122              :       detail::CompleteStep<FieldsTag, OptionsGroup, Label, ArraySectionIdTag>>;
     123              : };
     124              : 
     125              : }  // namespace newton_raphson
     126              : }  // namespace NonlinearSolver