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] function 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 : * @return a pair in which the first member is the value z such that f(x, z) = y and the second
39 : * member is dy/dz
40 : */
41 : template <typename T, typename Functor>
42 : std::pair<T, T>
43 197 : NewtonSolve(const T & x,
44 : const T & y,
45 : const Real z_initial_guess,
46 : const Real tolerance,
47 : const Functor & func,
48 : const std::string & caller_name,
49 : const unsigned int max_its = 100)
50 : {
51 : // R represents residual
52 :
53 : std::function<bool(const T &, const T &)> abs_tol_check =
54 : [tolerance](const T & R, const T & /*y*/)
55 9 : { return std::abs(MetaPhysicL::raw_value(R)) < tolerance; };
56 227 : std::function<bool(const T &, const T &)> rel_tol_check = [tolerance](const T & R, const T & y)
57 597 : { return std::abs(MetaPhysicL::raw_value(R / y)) < tolerance; };
58 393 : auto convergence_check = MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(y), 0, tolerance)
59 : ? abs_tol_check
60 : : rel_tol_check;
61 :
62 75 : T z = z_initial_guess, R, new_y, dy_dx, dy_dz;
63 : unsigned int iteration = 0;
64 :
65 : using std::isnan;
66 :
67 : do
68 : {
69 266 : func(x, z, new_y, dy_dx, dy_dz);
70 606 : R = new_y - y;
71 :
72 : // We always want to perform at least one update in order to get derivatives on z correct (z
73 : // corresponding to the initial guess will have no derivative information), so we don't
74 : // immediately return if we are converged
75 606 : const bool converged = convergence_check(R, y);
76 :
77 : #ifndef NDEBUG
78 : static constexpr Real perturbation_factor = 1 + 1e-8;
79 : T perturbed_y, dummy, dummy2;
80 : func(x, perturbation_factor * z, perturbed_y, dummy, dummy2);
81 : // Check the accuracy of the Jacobian
82 : auto J_differenced = (perturbed_y - new_y) / (1e-8 * z);
83 : if (!MooseUtils::relativeFuzzyEqual(J_differenced, dy_dz, 1e-2))
84 : mooseDoOnce(mooseWarning(caller_name + ": Bad Jacobian in NewtonSolve"));
85 : #endif
86 :
87 606 : z += -(R / dy_dz);
88 :
89 : // Check for NaNs
90 606 : if (isnan(z))
91 0 : mooseException(caller_name + ": NaN detected in Newton solve");
92 :
93 606 : if (converged)
94 : break;
95 409 : } while (++iteration < max_its);
96 :
97 : // Check for divergence or slow convergence of Newton's method
98 197 : if (iteration >= max_its)
99 0 : mooseException(caller_name +
100 : ": Newton solve convergence failed: maximum number of iterations, ",
101 : max_its,
102 : ", exceeded");
103 :
104 314 : return {z, dy_dz};
105 : }
106 :
107 : /**
108 : * NewtonSolve2D does a 2D Newton Solve to solve for the x and y such that:
109 : * f = func1(x, y) and g = func2(x, y). This is done for example in the constant of (v, e)
110 : * to (p, T) variable set conversion.
111 : * @param[in] f target value for func1
112 : * @param[in] g target value for func2
113 : * @param[in] x0 initial guess for first output variable
114 : * @param[in] y0 initial guess for second output variable
115 : * @param[out] x_final output for first variable
116 : * @param[out] y_final output for second variable
117 : * @param[in] f_tol criterion for relative or absolute (if f is sufficently close to zero)
118 : * convergence checking
119 : * @param[in] g_tol criterion for relative or absolute (if g is sufficently close to zero)
120 : * convergence checking
121 : * @param[in] func1 two-variable function returning both values and derivatives as references
122 : * @param[in] func2 two-variable function returning both values and derivatives as references
123 : * @param[in] max_its the maximum number of iterations for Newton's method
124 : */
125 : template <typename T, typename Functor1, typename Functor2>
126 : void
127 762478 : NewtonSolve2D(const T & f,
128 : const T & g,
129 : const Real x0,
130 : const Real y0,
131 : T & x_final,
132 : T & y_final,
133 : const Real f_tol,
134 : const Real g_tol,
135 : const Functor1 & func1,
136 : const Functor2 & func2,
137 : const unsigned int max_its = 100)
138 : {
139 :
140 : constexpr unsigned int system_size = 2;
141 762478 : DenseVector<T> targets = {{f, g}};
142 937001 : DenseVector<Real> tolerances = {{f_tol, g_tol}};
143 : // R represents a residual equal to y - y_in
144 20643455 : auto convergence_check = [&targets, &tolerances](const auto & minus_R)
145 : {
146 : using std::abs;
147 :
148 21413505 : for (const auto i : index_range(minus_R))
149 : {
150 20825571 : const auto error = abs(MooseUtils::absoluteFuzzyEqual(targets(i), 0, tolerances(i))
151 22 : ? minus_R(i)
152 20825527 : : minus_R(i) / targets(i));
153 20825549 : if (error >= tolerances(i))
154 : return false;
155 : }
156 : return true;
157 : };
158 :
159 937001 : DenseVector<T> u = {{x0, y0}};
160 762478 : DenseVector<T> minus_R(system_size), func_evals(system_size), u_update;
161 762478 : DenseMatrix<T> J(system_size, system_size);
162 : unsigned int iteration = 0;
163 : #ifndef NDEBUG
164 : DenseVector<Real> svs(system_size), evs_real(system_size), evs_imag(system_size);
165 : DenseMatrix<Real> raw_J(system_size, system_size), raw_J2(system_size, system_size);
166 : #endif
167 :
168 : typedef std::function<void(const T &, const T &, T &, T &, T &)> FuncType;
169 : std::array<FuncType, 2> func = {{func1, func2}};
170 :
171 12 : auto assign_solution = [&u, &x_final, &y_final]()
172 : {
173 762478 : x_final = u(0);
174 762478 : y_final = u(1);
175 : };
176 :
177 : using std::isnan, std::max, std::abs;
178 :
179 : do
180 : {
181 59642931 : for (const auto i : make_range(system_size))
182 39761954 : func[i](u(0), u(1), func_evals(i), J(i, 0), J(i, 1));
183 :
184 59642931 : for (const auto i : make_range(system_size))
185 39761954 : minus_R(i) = targets(i) - func_evals(i);
186 :
187 : // We always want to perform at least one update in order to get derivatives on z correct (z
188 : // corresponding to the initial guess will have no derivative information), so we don't
189 : // immediately return if we are converged
190 19880977 : const bool converged = convergence_check(minus_R);
191 :
192 : // Check for NaNs before proceeding to system solve. We may simultaneously not have NaNs in z
193 : // but have NaNs in the function evaluation
194 59642882 : for (const auto i : make_range(system_size))
195 39761930 : if (isnan(minus_R(i)))
196 : {
197 0 : assign_solution();
198 25 : mooseException("NaN detected in Newton solve");
199 : }
200 :
201 : // Do some Jacobi (rowmax) preconditioning
202 59642856 : for (const auto i : make_range(system_size))
203 : {
204 39761904 : const auto rowmax = max(abs(J(i, 0)), abs(J(i, 1)));
205 119285712 : for (const auto j : make_range(system_size))
206 79523808 : J(i, j) /= rowmax;
207 39761904 : minus_R(i) /= rowmax;
208 : }
209 :
210 : #ifndef NDEBUG
211 : //
212 : // Check nature of linearized system
213 : //
214 : for (const auto i : make_range(system_size))
215 : for (const auto j : make_range(system_size))
216 : {
217 : raw_J(i, j) = MetaPhysicL::raw_value(J(i, j));
218 : raw_J2(i, j) = MetaPhysicL::raw_value(J(i, j));
219 : }
220 : raw_J.svd(svs);
221 : raw_J2.evd(evs_real, evs_imag);
222 : #endif
223 :
224 19880952 : J.lu_solve(minus_R, u_update);
225 : // reset the decomposition
226 15 : J.zero();
227 19880952 : u += u_update;
228 :
229 : // Check for NaNs
230 59642664 : for (const auto i : make_range(system_size))
231 39761808 : if (isnan(u(i)))
232 : {
233 0 : assign_solution();
234 96 : mooseException("NaN detected in Newton solve");
235 : }
236 :
237 19880856 : if (converged)
238 : break;
239 19292901 : } while (++iteration < max_its);
240 :
241 4 : assign_solution();
242 :
243 : // Check for divergence or slow convergence of Newton's method
244 762357 : if (iteration >= max_its)
245 174402 : mooseException(
246 : "Newton solve convergence failed: maximum number of iterations, ", max_its, ", exceeded");
247 1350437 : }
248 : } // namespace FluidPropertiesUtils
|
