Line data Source code
1 : // Copyright 2024, UChicago Argonne, LLC
2 : // All Rights Reserved
3 : // Software Name: NEML2 -- the New Engineering material Model Library, version 2
4 : // By: Argonne National Laboratory
5 : // OPEN SOURCE LICENSE (MIT)
6 : //
7 : // Permission is hereby granted, free of charge, to any person obtaining a copy
8 : // of this software and associated documentation files (the "Software"), to deal
9 : // in the Software without restriction, including without limitation the rights
10 : // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11 : // copies of the Software, and to permit persons to whom the Software is
12 : // furnished to do so, subject to the following conditions:
13 : //
14 : // The above copyright notice and this permission notice shall be included in
15 : // all copies or substantial portions of the Software.
16 : //
17 : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18 : // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19 : // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20 : // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21 : // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22 : // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
23 : // THE SOFTWARE.
24 :
25 : #include <iostream>
26 : #include <iomanip>
27 :
28 : #include "neml2/solvers/NewtonWithLineSearch.h"
29 : #include "neml2/tensors/functions/bvv.h"
30 : #include "neml2/tensors/functions/sqrt.h"
31 : #include "neml2/tensors/assertions.h"
32 :
33 : namespace neml2
34 : {
35 : register_NEML2_object(NewtonWithLineSearch);
36 :
37 : OptionSet
38 5 : NewtonWithLineSearch::expected_options()
39 : {
40 5 : OptionSet options = Newton::expected_options();
41 5 : options.doc() = "The Newton-Raphson solver with line search.";
42 :
43 20 : EnumSelection linesearch_type({"BACKTRACKING", "STRONG_WOLFE"}, "BACKTRACKING");
44 5 : options.set<EnumSelection>("linesearch_type") = linesearch_type;
45 10 : options.set("linesearch_type").doc() = "The type of linesearch used."
46 10 : "Default: BACKTRACKING. Options are " +
47 15 : linesearch_type.candidates_str();
48 :
49 10 : options.set<unsigned int>("max_linesearch_iterations") = 10;
50 15 : options.set("max_linesearch_iterations").doc() =
51 : "Maximum allowable linesearch iterations. No error is produced upon reaching the maximum "
52 5 : "number of iterations, and the scale factor in the last iteration is used to scale the step.";
53 :
54 10 : options.set<double>("linesearch_cutback") = 2.0;
55 15 : options.set("linesearch_cutback").doc() = "Linesearch cut-back factor when the current scale "
56 5 : "factor cannot sufficiently reduce the residual.";
57 :
58 10 : options.set<double>("linesearch_stopping_criteria") = 1.0e-3;
59 15 : options.set("linesearch_stopping_criteria").doc() =
60 5 : "The lineseach tolerance slightly relaxing the definition of residual decrease";
61 :
62 10 : options.set<bool>("check_negative_critertia_value") = false;
63 10 : options.set("check_negative_critertia_value").doc() =
64 5 : "Whether to check if the convergence criteria for line search becomes negative";
65 :
66 10 : return options;
67 5 : }
68 :
69 4 : NewtonWithLineSearch::NewtonWithLineSearch(const OptionSet & options)
70 : : Newton(options),
71 8 : _linesearch_miter(options.get<unsigned int>("max_linesearch_iterations")),
72 8 : _linesearch_sigma(options.get<double>("linesearch_cutback")),
73 8 : _linesearch_c(options.get<double>("linesearch_stopping_criteria")),
74 8 : _type(options.get<EnumSelection>("linesearch_type")),
75 8 : _check_crit(options.get<bool>("check_negative_critertia_value"))
76 : {
77 4 : }
78 :
79 : void
80 79 : NewtonWithLineSearch::update(NonlinearSystem & system,
81 : NonlinearSystem::Sol<true> & x,
82 : const NonlinearSystem::Res<true> & r,
83 : const NonlinearSystem::Jac<true> & J)
84 : {
85 79 : auto dx = solve_direction(r, J);
86 79 : auto alpha = linesearch(system, x, dx, r);
87 79 : x = NonlinearSystem::Sol<true>(x.variable_data() + alpha * Tensor(dx));
88 79 : }
89 :
90 : Scalar
91 79 : NewtonWithLineSearch::linesearch(NonlinearSystem & system,
92 : const NonlinearSystem::Sol<true> & x,
93 : const NonlinearSystem::Sol<true> & dx,
94 : const NonlinearSystem::Res<true> & R0) const
95 : {
96 79 : auto alpha = Scalar::ones(x.batch_sizes(), x.options());
97 79 : const auto nR02 = bvv(R0, R0);
98 79 : auto crit = nR02;
99 :
100 192 : for (std::size_t i = 1; i < _linesearch_miter; i++)
101 : {
102 173 : NonlinearSystem::Sol<true> xp(Tensor(x) + alpha * Tensor(dx));
103 173 : auto R = system.residual(xp);
104 173 : auto nR2 = bvv(R, R);
105 :
106 346 : if (_type == "BACKTRACKING")
107 173 : crit = nR02 + 2.0 * _linesearch_c * alpha * bvv(R0, dx);
108 0 : else if (_type == "STRONG_WOLFE")
109 0 : crit = (1.0 - _linesearch_c * alpha) * nR02;
110 :
111 173 : if (verbose)
112 0 : std::cout << " LS ITERATION " << std::setw(3) << i << ", min(alpha) = " << std::scientific
113 0 : << at::min(alpha).item<double>() << ", max(||R||) = " << std::scientific
114 0 : << at::max(sqrt(nR2)).item<double>() << ", min(||Rc||) = " << std::scientific
115 0 : << at::min(sqrt(crit)).item<double>() << std::endl;
116 :
117 173 : auto stop = at::logical_or(nR2 <= crit, nR2 <= std::pow(atol, 2));
118 :
119 173 : if (at::all(stop).item<bool>())
120 60 : break;
121 :
122 226 : alpha.batch_index_put_({at::logical_not(stop)},
123 452 : alpha.batch_index({at::logical_not(stop)}) / _linesearch_sigma);
124 353 : }
125 :
126 79 : if (_check_crit)
127 0 : if (at::max(crit).item<double>() < 0)
128 : std::cerr << "WARNING: Line Search produces negative stopping "
129 : "criteria, this could lead to convergence issue. Try with other "
130 : "linesearch_type, increase linesearch_cutback "
131 0 : "or reduce linesearch_stopping_criteria"
132 0 : << std::endl;
133 :
134 158 : return alpha;
135 305 : }
136 :
137 : } // namespace neml2
|
