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
|
