NewtonRaphson.hpp
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 ObserveActions = tmpl::list<>,
97  typename Label = OptionsGroup>
98  using solve = tmpl::list<
99  ApplyNonlinearOperator,
100  detail::PrepareSolve<FieldsTag, OptionsGroup, Label>, ObserveActions,
102  detail::PrepareStep<FieldsTag, OptionsGroup, Label>, SolveLinearization,
103  detail::PerformStep<FieldsTag, OptionsGroup, Label>,
104  ApplyNonlinearOperator,
105  detail::ContributeToResidualMagnitudeReduction<FieldsTag, OptionsGroup,
106  Label>,
107  detail::Globalize<FieldsTag, OptionsGroup, Label>, ObserveActions,
108  detail::CompleteStep<FieldsTag, OptionsGroup, Label>>;
109 };
110
111 } // namespace newton_raphson
112 } // namespace NonlinearSolver
Tags.hpp