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 :
27 : #include "neml2/solvers/NonlinearSystem.h"
28 : #include "neml2/misc/assertions.h"
29 : #include "neml2/tensors/functions/bmm.h"
30 : #include "neml2/tensors/functions/diag_embed.h"
31 :
32 : namespace neml2
33 : {
34 : bool &
35 9184 : currently_solving_nonlinear_system()
36 : {
37 : thread_local bool _solving_nl_sys = false;
38 9184 : return _solving_nl_sys;
39 : }
40 :
41 199 : SolvingNonlinearSystem::SolvingNonlinearSystem(bool solving)
42 199 : : prev_bool(currently_solving_nonlinear_system())
43 : {
44 199 : currently_solving_nonlinear_system() = solving;
45 199 : }
46 :
47 199 : SolvingNonlinearSystem::~SolvingNonlinearSystem()
48 : {
49 199 : currently_solving_nonlinear_system() = prev_bool;
50 199 : }
51 :
52 : OptionSet
53 494 : NonlinearSystem::expected_options()
54 : {
55 494 : OptionSet options;
56 :
57 988 : options.set<bool>("automatic_scaling") = false;
58 1482 : options.set("automatic_scaling").doc() =
59 : "Whether to perform automatic scaling. See neml2::NonlinearSystem::init_scaling for "
60 494 : "implementation details.";
61 :
62 988 : options.set<double>("automatic_scaling_tol") = 0.01;
63 1482 : options.set("automatic_scaling_tol").doc() =
64 494 : "Tolerance used in iteratively updating the scaling matrices.";
65 :
66 988 : options.set<unsigned int>("automatic_scaling_miter") = 20;
67 988 : options.set("automatic_scaling_miter").doc() =
68 : "Maximum number of automatic scaling iterations. No error is produced upon reaching the "
69 : "maximum number of scaling iterations, and the scaling matrices obtained at the last "
70 494 : "iteration are used to scale the nonlinear system.";
71 :
72 494 : return options;
73 0 : }
74 :
75 : void
76 481 : NonlinearSystem::disable_automatic_scaling(OptionSet & options)
77 : {
78 962 : options.set("automatic_scaling").suppressed() = true;
79 962 : options.set("automatic_scaling_tol").suppressed() = true;
80 481 : options.set("automatic_scaling_miter").suppressed() = true;
81 481 : }
82 :
83 : void
84 12 : NonlinearSystem::enable_automatic_scaling(OptionSet & options)
85 : {
86 24 : options.set("automatic_scaling").suppressed() = false;
87 24 : options.set("automatic_scaling_tol").suppressed() = false;
88 12 : options.set("automatic_scaling_miter").suppressed() = false;
89 12 : }
90 :
91 494 : NonlinearSystem::NonlinearSystem(const OptionSet & options)
92 988 : : _autoscale(options.get<bool>("automatic_scaling")),
93 988 : _autoscale_tol(options.get<double>("automatic_scaling_tol")),
94 1482 : _autoscale_miter(options.get<unsigned int>("automatic_scaling_miter"))
95 : {
96 494 : }
97 :
98 : void
99 195 : NonlinearSystem::init_scaling(const NonlinearSystem::Sol<false> & x, const bool verbose)
100 : {
101 195 : if (!_autoscale)
102 191 : return;
103 :
104 4 : if (_scaling_matrices_initialized)
105 0 : return;
106 :
107 : // First compute the unscaled Jacobian
108 4 : auto J = Jacobian(x);
109 :
110 : // Initialize the scaling matrices
111 4 : auto Jp = J.clone();
112 4 : _row_scaling = Tensor::ones_like(x);
113 4 : _col_scaling = Tensor::ones_like(x);
114 :
115 4 : if (verbose)
116 0 : std::cout << "Before automatic scaling cond(J) = " << std::scientific
117 0 : << at::max(at::linalg_cond(Jp)).item<double>() << std::endl;
118 :
119 4 : const auto eps = machine_precision(x.scalar_type());
120 :
121 8 : for (unsigned int itr = 0; itr < _autoscale_miter; itr++)
122 : {
123 : // check for convergence
124 16 : auto rR = at::max(at::abs(1.0 - 1.0 / (at::sqrt(at::abs(std::get<0>(Jp.max(-1)))) + eps)))
125 8 : .item<double>();
126 16 : auto rC = at::max(at::abs(1.0 - 1.0 / (at::sqrt(at::abs(std::get<0>(Jp.max(-2)))) + eps)))
127 8 : .item<double>();
128 :
129 8 : if (verbose)
130 0 : std::cout << "ITERATION " << itr << ", ROW ILLNESS = " << std::scientific << rR
131 0 : << ", COL ILLNESS = " << std::scientific << rC << std::endl;
132 8 : if (rR < _autoscale_tol && rC < _autoscale_tol)
133 4 : break;
134 :
135 : // scale rows and columns
136 20 : for (Size i = 0; i < x.base_size(-1); i++)
137 : {
138 32 : auto ar = 1.0 / (at::sqrt(at::max(at::abs(Jp.base_index({i})))) + eps);
139 64 : auto ac = 1.0 / (at::sqrt(at::max(at::abs(Jp.base_index({indexing::Slice(), i})))) + eps);
140 32 : _row_scaling.base_index({i}) *= ar;
141 32 : _col_scaling.base_index({i}) *= ac;
142 32 : Jp.base_index({i}) *= ar;
143 64 : Jp.base_index({indexing::Slice(), i}) *= ac;
144 16 : }
145 : }
146 :
147 4 : _scaling_matrices_initialized = true;
148 :
149 4 : if (verbose)
150 0 : std::cout << " After automatic scaling cond(J) = " << std::scientific
151 0 : << at::max(at::linalg_cond(Jp)).item<double>() << std::endl;
152 132 : }
153 :
154 : void
155 98 : NonlinearSystem::ensure_scaling_matrices_initialized_dbg() const
156 : {
157 98 : neml_assert_dbg(
158 98 : _autoscale == _scaling_matrices_initialized,
159 98 : _autoscale ? "Automatic scaling is requested but scaling matrices have not been initialized."
160 : : "Automatic scaling is not requested but scaling matrices were initialized.");
161 98 : }
162 :
163 : NonlinearSystem::Res<true>
164 1737 : NonlinearSystem::scale(const NonlinearSystem::Res<false> & r) const
165 : {
166 1737 : if (!_autoscale)
167 1702 : return Res<true>(r);
168 :
169 35 : ensure_scaling_matrices_initialized_dbg();
170 35 : return Res<true>(_row_scaling * r);
171 : }
172 :
173 : NonlinearSystem::Res<false>
174 0 : NonlinearSystem::unscale(const NonlinearSystem::Res<true> & r) const
175 : {
176 0 : if (!_autoscale)
177 0 : return Res<false>(r);
178 :
179 0 : ensure_scaling_matrices_initialized_dbg();
180 0 : return Res<false>(1. / _row_scaling * r);
181 : }
182 :
183 : NonlinearSystem::Jac<true>
184 1355 : NonlinearSystem::scale(const NonlinearSystem::Jac<false> & J) const
185 : {
186 1355 : if (!_autoscale)
187 1335 : return Jac<true>(J);
188 :
189 20 : ensure_scaling_matrices_initialized_dbg();
190 20 : return Jac<true>(bmm(bmm(base_diag_embed(_row_scaling), J), base_diag_embed(_col_scaling)));
191 : }
192 :
193 : NonlinearSystem::Jac<false>
194 0 : NonlinearSystem::unscale(const NonlinearSystem::Jac<true> & J) const
195 : {
196 0 : if (!_autoscale)
197 0 : return Jac<false>(J);
198 :
199 0 : ensure_scaling_matrices_initialized_dbg();
200 : return Jac<false>(
201 0 : bmm(bmm(base_diag_embed(1.0 / _row_scaling), J), base_diag_embed(1.0 / _col_scaling)));
202 : }
203 :
204 : NonlinearSystem::Sol<true>
205 219 : NonlinearSystem::scale(const NonlinearSystem::Sol<false> & u) const
206 : {
207 219 : if (!_autoscale)
208 215 : return Sol<true>(u);
209 :
210 4 : ensure_scaling_matrices_initialized_dbg();
211 4 : return Sol<true>(1. / _col_scaling * u);
212 : }
213 :
214 : NonlinearSystem::Sol<false>
215 1956 : NonlinearSystem::unscale(const NonlinearSystem::Sol<true> & u) const
216 : {
217 1956 : if (!_autoscale)
218 1917 : return Sol<false>(u);
219 :
220 39 : ensure_scaling_matrices_initialized_dbg();
221 39 : return Sol<false>(_col_scaling * u);
222 : }
223 :
224 : void
225 1738 : NonlinearSystem::set_guess(const NonlinearSystem::Sol<true> & x)
226 : {
227 1738 : set_guess(unscale(x));
228 1738 : }
229 : } // namespace neml2
|
