Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : // * This file is part of the MOOSE framework
11 : // * https://mooseframework.inl.gov
12 : // *
13 : // * All rights reserved, see COPYRIGHT for full restrictions
14 : // * https://github.com/idaholab/moose/blob/master/COPYRIGHT
15 : // *
16 : // * Licensed under LGPL 2.1, please see LICENSE for details
17 : // * https://www.gnu.org/licenses/lgpl-2.1.html
18 :
19 : #pragma once
20 :
21 : #include "Moose.h"
22 : #include "MooseUtils.h"
23 : #include "libmesh/dense_matrix.h"
24 :
25 : namespace FluidPropertiesUtils
26 : {
27 : /**
28 : * NewtonSolve does a 1D Newton Solve to solve the equation y = f(x, z) for variable z.
29 : * @param[in] x constant first argument of the f(x, z) term
30 : * @param[in] y constant which should be equal to f(x, z) with a converged z
31 : * @param[in] z_initial_guess initial guess for return variables
32 : * @param[in] tolerance criterion for relative or absolute (if y is sufficiently close to zero)
33 : * convergence checking
34 : * @param[in] y_from_x_z two-variable function returning both values and derivatives as references
35 : * @param[in] caller_name name of the fluid properties appended to name of the routine calling the
36 : * method
37 : * @param[in] max_its the maximum number of iterations for Newton's method
38 : * @param[in] verbose whether to output Newton iteration data
39 : * @return a pair in which the first member is the value z such that f(x, z) = y and the second
40 : * member is dy/dz
41 : */
42 : template <typename T, typename Functor>
43 : std::pair<T, T>
44 456 : NewtonSolve(const T & x,
45 : const T & y,
46 : const Real z_initial_guess,
47 : const Real tolerance,
48 : const Functor & y_from_x_z,
49 : const std::string & caller_name,
50 : const unsigned int max_its = 100,
51 : const bool verbose = false)
52 : {
53 : // R represents residual
54 :
55 : std::function<bool(const T &, const T &)> abs_tol_check =
56 : [tolerance](const T & R, const T & /*y*/)
57 9 : { return std::abs(MetaPhysicL::raw_value(R)) < tolerance; };
58 295 : std::function<bool(const T &, const T &)> rel_tol_check = [tolerance](const T & R, const T & y)
59 1104 : { return std::abs(MetaPhysicL::raw_value(R / y)) < tolerance; };
60 911 : auto convergence_check = MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(y), 0, tolerance)
61 : ? abs_tol_check
62 : : rel_tol_check;
63 :
64 93 : T z = z_initial_guess, R, new_y, dy_dx, dy_dz;
65 : unsigned int iteration = 0;
66 :
67 : using std::isnan;
68 456 : if (verbose)
69 : Moose::out << "Target value for 1D Newton inversion:\n" << y << std::endl;
70 :
71 : do
72 : {
73 693 : y_from_x_z(x, z, new_y, dy_dx, dy_dz);
74 1113 : R = new_y - y;
75 :
76 : // We always want to perform at least one update in order to get derivatives on z correct (z
77 : // corresponding to the initial guess will have no derivative information), so we don't
78 : // immediately return if we are converged
79 1113 : const bool converged = convergence_check(R, y);
80 :
81 : #ifndef NDEBUG
82 : static constexpr Real perturbation_factor = 1 + 1e-8;
83 : T perturbed_y, dummy, dummy2;
84 : y_from_x_z(x, perturbation_factor * z, perturbed_y, dummy, dummy2);
85 : // Check the accuracy of the Jacobian
86 : auto J_differenced = (perturbed_y - new_y) / (1e-8 * z);
87 : if (!MooseUtils::relativeFuzzyEqual(J_differenced, dy_dz, 1e-2))
88 : mooseDoOnce(mooseWarning(caller_name + ": Bad Jacobian in NewtonSolve"));
89 : #endif
90 :
91 1113 : z += -(R / dy_dz);
92 :
93 1113 : if (verbose)
94 : {
95 : Moose::out << "Iteration " << iteration << std::endl;
96 : Moose::out << "Current solution vector: " << z << std::endl;
97 0 : Moose::out << "Current (minus) residual: " << -R << std::endl;
98 : Moose::out << "Current Jacobian: " << dy_dz << std::endl;
99 : }
100 :
101 : // Check for NaNs
102 1113 : if (isnan(z))
103 0 : mooseException(caller_name + ": NaN detected in Newton solve");
104 :
105 1113 : if (converged)
106 : break;
107 657 : } while (++iteration < max_its);
108 :
109 : // Check for divergence or slow convergence of Newton's method
110 456 : if (iteration >= max_its)
111 0 : mooseException(caller_name +
112 : ": Newton solve convergence failed: maximum number of iterations, ",
113 : max_its,
114 : ", exceeded");
115 :
116 : // z was updated, we need to recompute the derivative
117 341 : y_from_x_z(x, z, new_y, dy_dx, dy_dz);
118 :
119 807 : return {z, dy_dz};
120 : }
121 :
122 : /**
123 : * NewtonSolve2D does a 2D Newton Solve to solve for the x and y such that:
124 : * f = f_from_x_y(x, y) and g = g_from_x_y(x, y). This is done for example in the constant of (v, e)
125 : * to (p, T) variable set conversion.
126 : * @param[in] f target value for f_from_x_y
127 : * @param[in] g target value for g_from_x_y
128 : * @param[in] x0 initial guess for first output variable
129 : * @param[in] y0 initial guess for second output variable
130 : * @param[out] x_final output for first variable
131 : * @param[out] y_final output for second variable
132 : * @param[in] f_tol criterion for relative or absolute (if f is sufficently close to zero)
133 : * convergence checking
134 : * @param[in] g_tol criterion for relative or absolute (if g is sufficently close to zero)
135 : * convergence checking
136 : * @param[in] f_from_x_y two-variable function returning both values and derivatives as references
137 : * @param[in] g_from_x_y two-variable function returning both values and derivatives as references
138 : * @param[in] caller_name routine calling this solve
139 : * @param[in] max_its the maximum number of iterations for Newton's method
140 : * @param[in] debug whether to output the solution, residual and Jacobian on every iteration
141 : */
142 : template <typename T, typename Functor1, typename Functor2>
143 : void
144 520395 : NewtonSolve2D(const T & f,
145 : const T & g,
146 : const Real x0,
147 : const Real y0,
148 : T & x_final,
149 : T & y_final,
150 : const Real f_tol,
151 : const Real g_tol,
152 : const Functor1 & f_from_x_y,
153 : const Functor2 & g_from_x_y,
154 : const std::string & caller_name = "",
155 : const unsigned int max_its = 100,
156 : bool debug = false)
157 : {
158 :
159 : constexpr unsigned int system_size = 2;
160 520395 : DenseVector<T> targets = {{f, g}};
161 640070 : DenseVector<Real> tolerances = {{f_tol, g_tol}};
162 : // R represents a residual equal to y - y_in
163 14138033 : auto convergence_check = [&targets, &tolerances](const auto & minus_R)
164 : {
165 : using std::abs;
166 :
167 14665736 : for (const auto i : index_range(minus_R))
168 : {
169 14265047 : const auto error = abs(MooseUtils::absoluteFuzzyEqual(targets(i), 0, tolerances(i))
170 32 : ? minus_R(i)
171 14264983 : : minus_R(i) / targets(i));
172 14265015 : if (error >= tolerances(i))
173 : return false;
174 : }
175 : return true;
176 : };
177 :
178 640070 : DenseVector<T> u = {{x0, y0}};
179 520395 : DenseVector<T> minus_R(system_size), func_evals(system_size), u_update(system_size);
180 520395 : DenseMatrix<T> J(system_size, system_size);
181 : unsigned int iteration = 0;
182 : #ifndef NDEBUG
183 : DenseVector<Real> svs(system_size), evs_real(system_size), evs_imag(system_size);
184 : DenseMatrix<Real> raw_J(system_size, system_size), raw_J2(system_size, system_size);
185 : #endif
186 :
187 : typedef std::function<void(const T &, const T &, T &, T &, T &)> FuncType;
188 : std::array<FuncType, 2> func = {{f_from_x_y, g_from_x_y}};
189 :
190 18 : auto assign_solution = [&u, &x_final, &y_final]()
191 : {
192 520395 : x_final = u(0);
193 520395 : y_final = u(1);
194 : };
195 759711 : auto status_string = [&u, &func_evals, &targets, &minus_R](unsigned int comp) -> std::stringstream
196 : {
197 239316 : std::stringstream ss;
198 478632 : ss << "Current solution for component " << comp << ": " << u(comp)
199 239316 : << " (current ordinate: " << func_evals(comp) << " -> target: " << targets(comp)
200 239316 : << ", scaled residual: " << minus_R(comp) << ")";
201 239316 : return ss;
202 0 : };
203 520395 : if (debug)
204 : Moose::out << "Target values for 2D Newton inversion:\n" << targets << std::endl;
205 :
206 : using std::isnan, std::max, std::abs;
207 :
208 : do
209 : {
210 40852914 : for (const auto i : make_range(system_size))
211 27235276 : func[i](u(0), u(1), func_evals(i), J(i, 0), J(i, 1));
212 :
213 40852914 : for (const auto i : make_range(system_size))
214 27235276 : minus_R(i) = targets(i) - func_evals(i);
215 :
216 : // We always want to perform at least one update in order to get derivatives on z correct (z
217 : // corresponding to the initial guess will have no derivative information), so we don't
218 : // immediately return if we are converged
219 13617638 : const bool converged = convergence_check(minus_R);
220 :
221 : // Check for NaNs before proceeding to system solve. We may simultaneously not have NaNs in z
222 : // but have NaNs in the function evaluation
223 40852881 : for (const auto i : make_range(system_size))
224 27235260 : if (isnan(minus_R(i)))
225 : {
226 0 : assign_solution();
227 34 : mooseException(caller_name + ": NaN detected in Newton solve");
228 : }
229 :
230 13617621 : if (debug)
231 : {
232 : Moose::out << "Iteration " << iteration << std::endl;
233 : Moose::out << "Current solution vector:\n" << u << std::endl;
234 : Moose::out << "Current (minus) residual:\n" << minus_R << std::endl;
235 : Moose::out << "Current Jacobian:\n" << J << std::endl;
236 : }
237 :
238 : // Do some Jacobi (rowmax) preconditioning and check for an empty row
239 : int degenerate_row = -1;
240 40852863 : for (const auto i : make_range(system_size))
241 : {
242 27235242 : const auto rowmax = max(abs(J(i, 0)), abs(J(i, 1)));
243 27235242 : if (rowmax > 0)
244 : {
245 81689520 : for (const auto j : make_range(system_size))
246 54459680 : J(i, j) /= rowmax;
247 27229840 : minus_R(i) /= rowmax;
248 : }
249 : else
250 : {
251 5402 : if (degenerate_row != -1)
252 0 : mooseException(caller_name + ": Jacobian is all zeros in NewtonSolve2D");
253 5402 : degenerate_row = i;
254 : }
255 : }
256 :
257 : #ifndef NDEBUG
258 : //
259 : // Check nature of linearized system
260 : //
261 : for (const auto i : make_range(system_size))
262 : for (const auto j : make_range(system_size))
263 : {
264 : raw_J(i, j) = MetaPhysicL::raw_value(J(i, j));
265 : raw_J2(i, j) = MetaPhysicL::raw_value(J(i, j));
266 : }
267 : raw_J.svd(svs);
268 : raw_J2.evd(evs_real, evs_imag);
269 : if (debug)
270 : Moose::out << "Jacobian singular values:\n" << svs << std::endl;
271 : #endif
272 :
273 13617621 : if (degenerate_row == -1)
274 13612219 : J.lu_solve(minus_R, u_update);
275 : else
276 : {
277 : // use a 1D newton when the Jacobian has an empty row
278 5402 : const auto other_row = system_size - 1 - degenerate_row;
279 5402 : u_update(other_row) = minus_R(other_row) / J(other_row, other_row);
280 5402 : u_update(degenerate_row) = 0;
281 : }
282 : // reset the decomposition
283 23 : J.zero();
284 13617621 : u += u_update;
285 :
286 : // Check for NaNs
287 40852863 : for (const auto i : make_range(system_size))
288 27235242 : if (isnan(u(i)))
289 : {
290 0 : assign_solution();
291 0 : mooseException(caller_name + ": NaN detected in NewtonSolve2D\n" + status_string(0).str() +
292 : "\n" + status_string(1).str());
293 : }
294 :
295 13617621 : if (converged)
296 : break;
297 13216901 : } while (++iteration < max_its);
298 :
299 6 : assign_solution();
300 :
301 : // Check for divergence or slow convergence of Newton's method
302 520378 : if (iteration >= max_its)
303 717948 : mooseException(caller_name +
304 : ": Newton solve convergence failed: maximum number of iterations, ",
305 : max_its,
306 : ", exceeded.\n" + status_string(0).str() + "\n" + status_string(1).str());
307 921121 : }
308 : } // namespace FluidPropertiesUtils
|
