NewtonRaphson.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include "DataStructures/DataBox/PrefixHelpers.hpp"
8 #include "IO/Observer/Helpers.hpp"
9 #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/ElementActions.hpp"
10 #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/ResidualMonitor.hpp"
11 #include "ParallelAlgorithms/NonlinearSolver/Observe.hpp"
13 #include "Utilities/TMPL.hpp"
14 
15 namespace NonlinearSolver {
16 /// Items related to the NewtonRaphson nonlinear solver
17 namespace newton_raphson {
18 
19 /*!
20  * \brief A Newton-Raphson correction scheme for nonlinear systems of equations
21  * \f$A_\mathrm{nonlinear}(x)=b\f$.
22  *
23  * We can use a correction scheme to solve a nonlinear problem
24  * \f$A_\mathrm{nonlinear}(x)=b\f$ by repeatedly solving a linearization of it.
25  * A Newton-Raphson scheme iteratively refines an initial guess \f$x_0\f$ by
26  * repeatedly solving the linearized problem
27  *
28  * \f{equation}
29  * \frac{\delta A_\mathrm{nonlinear}}{\delta x}(x_k) \delta x_k =
30  * b-A_\mathrm{nonlinear}(x_k) \equiv r_k
31  * \f}
32  *
33  * for the correction \f$\delta x_k\f$ and then updating the solution as
34  * \f$x_{k+1}=x_k + \delta x_k\f$.
35  *
36  * The operations we need to supply to the algorithm are the nonlinear operator
37  * \f$A_\mathrm{nonlinear}(x)\f$ and a linear solver for the linearized problem
38  * \f$\frac{\delta A_\mathrm{nonlinear}}{\delta x}(x_k) \delta x_k = r_k\f$.
39  * Each step of the algorithm expects that \f$A_\mathrm{nonlinear}(x)\f$ is
40  * computed and stored in the DataBox as
41  * `db::add_tag_prefix<NonlinearSolver::Tags::OperatorAppliedTo, FieldsTag>`.
42  * To perform a solve, add the `solve` action list to an array parallel
43  * component. Pass the actions that compute \f$A_\mathrm{nonlinear}(x)\f$ as
44  * the first template parameter to `solve`. As the second template parameter,
45  * pass the action list that performs a linear solve of the linearized operator
46  * \f$\frac{\delta A_\mathrm{nonlinear}}{\delta x}(x)\f$ for the field
47  * \f$\delta x\f$ (the `linear_solver_fields_tag`) sourced by
48  * \f$r\f$ (the `linear_solver_source_tag`). You will find suitable iterative
49  * linear solvers in the `LinearSolver` namespace. The third template parameter
50  * allows you to pass any further actions you wish to run in each step of the
51  * algorithm (such as observations). If you add the `solve` action list multiple
52  * times, use the fourth template parameter to label each solve with a different
53  * type.
54  *
55  * \par Globalization:
56  * This nonlinear solver supports a line-search (or "backtracking")
57  * globalization. If a step does not sufficiently decrease the residual (see
58  * `NonlinearSolver::OptionTags::SufficientDecrease` for details on the
59  * sufficient decrease condition), the step length is reduced until the residual
60  * sufficiently decreases or the maximum number of globalization steps is
61  * reached. The reduction in step length is determined by the minimum of a
62  * quadratic (first globalization step) or cubic (subsequent globalization
63  * steps) polynomial interpolation according to Algorithm A6.1.3 in
64  * \cite DennisSchnabel (p. 325) (see
65  * `NonlinearSolver::newton_raphson::next_step_length`). Alternatives to a
66  * line-search globalization, such as a trust-region globalization or more
67  * sophisticated nonlinear preconditioning techniques (see e.g. \cite Brune2015
68  * for an overview), are not currently implemented.
69  */
70 template <typename Metavariables, typename FieldsTag, typename OptionsGroup,
71  typename SourceTag =
73 struct NewtonRaphson {
74  using fields_tag = FieldsTag;
75  using options_group = OptionsGroup;
76  using source_tag = SourceTag;
77 
78  using operand_tag = fields_tag;
79  using linear_solver_fields_tag =
81  using linear_solver_source_tag =
83 
84  using component_list = tmpl::list<
85  detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
86 
87  using initialize_element =
88  detail::InitializeElement<FieldsTag, OptionsGroup, SourceTag>;
89 
90  using register_element = tmpl::list<>;
91 
92  using observed_reduction_data_tags = observers::make_reduction_data_tags<
93  tmpl::list<NonlinearSolver::observe_detail::reduction_data>>;
94 
95  template <typename ApplyNonlinearOperator, typename SolveLinearization,
96  typename CompleteStepActions = tmpl::list<>,
97  typename Label = OptionsGroup>
98  using solve = tmpl::list<
99  detail::PrepareSolve<FieldsTag, OptionsGroup, Label>,
100  detail::ReceiveInitialHasConverged<FieldsTag, OptionsGroup, Label>,
101  detail::PrepareStep<FieldsTag, OptionsGroup, Label>, SolveLinearization,
102  detail::PerformStep<FieldsTag, OptionsGroup, Label>,
103  ApplyNonlinearOperator,
104  detail::ContributeToResidualMagnitudeReduction<FieldsTag, OptionsGroup,
105  Label>,
106  detail::Globalize<FieldsTag, OptionsGroup, Label>, CompleteStepActions,
107  detail::CompleteStep<FieldsTag, OptionsGroup, Label>>;
108 };
109 
110 } // namespace newton_raphson
111 } // namespace NonlinearSolver
Tags.hpp
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
RootFinder::newton_raphson
DataVector newton_raphson(const Function &f, const DataVector &initial_guess, const DataVector &lower_bound, const DataVector &upper_bound, const size_t digits, const size_t max_iterations=50)
Finds the root of the function f with the Newton-Raphson method on each element in a DataVector.
Definition: NewtonRaphson.hpp:89
NonlinearSolver
Functionality for solving nonlinear systems of equations.
Definition: ElementActions.hpp:49
Prefixes.hpp
NonlinearSolver::newton_raphson::NewtonRaphson
A Newton-Raphson correction scheme for nonlinear systems of equations .
Definition: NewtonRaphson.hpp:73
TMPL.hpp