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