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 : #pragma once
11 :
12 : #include "MooseTypes.h"
13 : #include "RankTwoTensor.h"
14 : #include "InputParameters.h"
15 :
16 : #include "libmesh/utility.h"
17 : #include "libmesh/int_range.h"
18 : #include "libmesh/vector_value.h"
19 :
20 : #include "Eigen/Core"
21 : #include <Eigen/Dense>
22 : #include <unsupported/Eigen/NonLinearOptimization>
23 :
24 : #include "EigenADReal.h"
25 :
26 : /**
27 : * Internal use namespace that contains template helpers to map each residual
28 : * type to its corresponding Jacobian type.
29 : */
30 : namespace NestedSolveInternal
31 : {
32 : template <typename T>
33 : struct CorrespondingJacobianTempl;
34 : }
35 :
36 : /**
37 : * NestedSolveTempl<is_ad> and its instantiations NestedSolve and ADNestedSolve
38 : * are utility classes that implement a non-linear solver. These classes can be
39 : * instantiated in any MOOSE object to perform local Newton-Raphson solves.
40 : */
41 : template <bool is_ad>
42 : class NestedSolveTempl
43 : {
44 : public:
45 : static InputParameters validParams();
46 :
47 : NestedSolveTempl();
48 : NestedSolveTempl(const InputParameters & params);
49 :
50 : ///@{ AD/non-AD switched type shortcuts
51 : using NSReal = Moose::GenericType<Real, is_ad>;
52 : using NSRealVectorValue = Moose::GenericType<RealVectorValue, is_ad>;
53 : using NSRankTwoTensor = Moose::GenericType<RankTwoTensor, is_ad>;
54 : ///@}
55 :
56 : ///@{ Eigen type shortcuts
57 : using DynamicVector = Eigen::Matrix<NSReal, Eigen::Dynamic, 1>;
58 : using DynamicMatrix = Eigen::Matrix<NSReal, Eigen::Dynamic, Eigen::Dynamic>;
59 : ///@}
60 :
61 : ///@{ Deduce the Jacobian type from the solution type
62 : template <typename T>
63 : struct CorrespondingJacobianTempl;
64 : template <typename T>
65 : using CorrespondingJacobian = typename NestedSolveInternal::CorrespondingJacobianTempl<T>::type;
66 : ///@}
67 :
68 : /// Residual and solution type
69 : template <int N = 0>
70 : using Value =
71 : typename std::conditional<N == 1,
72 : NSReal,
73 : typename std::conditional<N == 0,
74 : NestedSolveTempl<is_ad>::DynamicVector,
75 : Eigen::Matrix<NSReal, N, 1>>::type>::type;
76 :
77 : /// Jacobian matrix type
78 : template <int N = 0>
79 : using Jacobian =
80 : typename std::conditional<N == 1,
81 : NSReal,
82 : typename std::conditional<N == 0,
83 : NestedSolveTempl<is_ad>::DynamicMatrix,
84 : Eigen::Matrix<NSReal, N, N>>::type>::type;
85 :
86 : /// Solve the N*N nonlinear equation system using a built-in Netwon-Raphson loop
87 : template <typename V, typename T>
88 : void nonlinear(V & guess, T && compute);
89 :
90 : /// Solve the N*N nonlinear equation system using the damped Netwon-Raphson loop
91 : template <typename V, typename T, typename C>
92 : void nonlinearDamped(V & guess, T && compute, C && computeCondition);
93 :
94 : ///@{ The separate residual/Jacobian functor versions use Eigen::HybridNonLinearSolver
95 : /// with a custom backwards compatible convergence check that allows for looser tolerances.
96 : template <typename R, typename J>
97 : void nonlinear(DynamicVector & guess, R & computeResidual, J & computeJacobian);
98 : template <typename R, typename J>
99 : void nonlinear(NSReal & guess, R & computeResidual, J & computeJacobian);
100 : template <typename R, typename J>
101 : void nonlinear(NSRealVectorValue & guess, R & computeResidual, J & computeJacobian);
102 : ///@}
103 :
104 : ///@{ The separate residual/Jacobian functor versions use Eigen::HybridNonLinearSolver
105 : /// with the built-in convergence check to tight numerical tolerance.
106 : template <typename R, typename J>
107 : void nonlinearTight(DynamicVector & guess, R & computeResidual, J & computeJacobian);
108 : template <typename R, typename J>
109 : void nonlinearTight(NSReal & guess, R & computeResidual, J & computeJacobian);
110 : template <typename R, typename J>
111 : void nonlinearTight(NSRealVectorValue & guess, R & computeResidual, J & computeJacobian);
112 : ///@}
113 :
114 : ///@{ Perform a bounded solve use Eigen::HybridNonLinearSolver
115 : template <typename R, typename J, typename B>
116 : void
117 : nonlinearBounded(NSReal & guess, R & computeResidual, J & computeJacobian, B & computeBounds);
118 : ///@}
119 : ///@{ default values
120 14575 : static Real relativeToleranceDefault() { return 1e-8; }
121 14575 : static Real absoluteToleranceDefault() { return 1e-13; }
122 14575 : static Real xToleranceDefault() { return 1e-15; }
123 14575 : static unsigned int minIterationsDefault() { return 3; }
124 14575 : static unsigned int maxIterationsDefault() { return 1000; }
125 14575 : static Real acceptableMultiplierDefault() { return 10.0; }
126 14575 : static Real dampingFactorDefault() { return 0.8; }
127 14575 : static unsigned int maxDampingIterationsDefault() { return 100; }
128 : ///@}
129 :
130 3 : void setRelativeTolerance(Real rel) { _relative_tolerance_square = rel * rel; }
131 0 : void setAbsoluteTolerance(Real abs) { _absolute_tolerance_square = abs * abs; }
132 0 : void setAcceptableMultiplier(Real am) { _acceptable_multiplier = am; }
133 :
134 : Real _relative_tolerance_square;
135 : Real _absolute_tolerance_square;
136 : /// Threshold for minimum step size of linear iterations
137 : Real _delta_thresh;
138 : /// Damping factor
139 : Real _damping_factor;
140 : unsigned int _max_damping_iterations;
141 :
142 : unsigned int _min_iterations;
143 : unsigned int _max_iterations;
144 : Real _acceptable_multiplier;
145 :
146 : /// possible solver states
147 : enum class State
148 : {
149 : NONE,
150 : CONVERGED_ABS,
151 : CONVERGED_REL,
152 : CONVERGED_ACCEPTABLE_ABS,
153 : CONVERGED_ACCEPTABLE_REL,
154 : CONVERGED_BOUNDS,
155 : EXACT_GUESS,
156 : CONVERGED_XTOL,
157 : NOT_CONVERGED
158 : };
159 :
160 : /// Get the solver state
161 1 : const State & getState() const { return _state; }
162 : /// Get the number of iterations from the last solve
163 0 : const std::size_t & getIterations() { return _n_iterations; };
164 :
165 : ///@{ Compute squared norm of v (dropping derivatives as this is only for convergence checking)
166 : template <typename V>
167 : static Real normSquare(const V & v);
168 : static Real normSquare(const NSReal & v);
169 : static Real normSquare(const NSRealVectorValue & v);
170 : ///@}
171 :
172 : ///@{ Check if |a| < |b * c| for all elements in a and b. This checks if 'a' is small relative to
173 : //'b' by a factor of 'c'
174 : template <typename V>
175 : static bool isRelSmall(const V & a, const V & b, const V & c);
176 : static bool isRelSmall(const NSReal & a, const NSReal & b, const NSReal & c);
177 : static bool
178 : isRelSmall(const NSRealVectorValue & a, const NSRealVectorValue & b, const NSReal & c);
179 : static bool isRelSmall(const DynamicVector & a, const DynamicVector & b, const NSReal & c);
180 : ///@}
181 :
182 : protected:
183 : /// current solver state
184 : State _state;
185 :
186 : /// number of nested iterations
187 : std::size_t _n_iterations;
188 :
189 : /// Size a dynamic Jacobian matrix correctly
190 : void sizeItems(const NestedSolveTempl<is_ad>::DynamicVector & guess,
191 : NestedSolveTempl<is_ad>::DynamicVector & residual,
192 : NestedSolveTempl<is_ad>::DynamicMatrix & jacobian) const;
193 :
194 : /// Sizing is a no-op for compile time sized types (and scalars)
195 : template <typename V, typename T>
196 7625 : void sizeItems(const V &, V &, T &) const
197 : {
198 7625 : }
199 :
200 : ///@{ Solve A*x=b for x
201 : template <typename J, typename V>
202 : void linear(const J & A, V & x, const V & b) const;
203 113626 : void linear(const NSReal & A, NSReal & x, const NSReal & b) const { x = b / A; }
204 : void linear(const NSRankTwoTensor & A, NSRealVectorValue & x, const NSRealVectorValue & b) const;
205 : ///@}
206 :
207 : /// Convergence check (updates _status)
208 : bool isConverged(Real r0_square, Real r_square, bool acceptable);
209 :
210 : private:
211 : ///@{ Build a suitable Eigen adaptor functors to interface between moose and Eigen types
212 : template <bool store_residual_norm, typename R, typename J>
213 : auto make_adaptor(R & residual, J & jacobian);
214 : template <bool store_residual_norm, typename R, typename J>
215 : auto make_real_adaptor(R & residual, J & jacobian);
216 : template <bool store_residual_norm, typename R, typename J>
217 : auto make_realvector_adaptor(R & residual, J & jacobian);
218 : ///@}
219 : };
220 :
221 : template <bool is_ad>
222 : template <typename R, typename J>
223 : void
224 1 : NestedSolveTempl<is_ad>::nonlinear(NestedSolveTempl<is_ad>::DynamicVector & guess,
225 : R & computeResidual,
226 : J & computeJacobian)
227 : {
228 : // adaptor functor to utilize the Eigen non-linear solver
229 1 : auto functor = make_adaptor<true>(computeResidual, computeJacobian);
230 1 : Eigen::HybridNonLinearSolver<decltype(functor), NSReal> solver(functor);
231 1 : auto status = solver.solve(guess);
232 :
233 1 : if (status == Eigen::HybridNonLinearSolverSpace::RelativeErrorTooSmall)
234 1 : _state = State::CONVERGED_REL;
235 : else
236 0 : _state = State::NOT_CONVERGED;
237 1 : }
238 :
239 : template <bool is_ad>
240 : template <typename R, typename J>
241 : void
242 1 : NestedSolveTempl<is_ad>::nonlinear(NSReal & guess, R & computeResidual, J & computeJacobian)
243 : {
244 : using V = NestedSolveTempl<is_ad>::DynamicVector;
245 :
246 : // adaptor functor to utilize the Eigen non-linear solver
247 1 : auto functor = make_real_adaptor<true>(computeResidual, computeJacobian);
248 1 : Eigen::HybridNonLinearSolver<decltype(functor), NSReal> solver(functor);
249 1 : V guess_eigen(1);
250 1 : guess_eigen << guess;
251 1 : const auto & r_square = functor.getResidualNorm();
252 1 : solver.solveInit(guess_eigen);
253 :
254 : // compute first residual norm for relative convergence checks
255 1 : auto r0_square = r_square;
256 1 : if (r0_square == 0)
257 : {
258 0 : _state = State::EXACT_GUESS;
259 0 : return;
260 : }
261 :
262 : // perform non-linear iterations
263 1 : _n_iterations = 1;
264 2 : while (_n_iterations <= _max_iterations &&
265 1 : solver.solveOneStep(guess_eigen) == Eigen::HybridNonLinearSolverSpace::Running)
266 : {
267 : // check convergence
268 0 : if (_n_iterations >= _min_iterations && isConverged(r0_square, r_square, /*acceptable=*/false))
269 : {
270 0 : guess = guess_eigen(0);
271 0 : return;
272 : }
273 0 : _n_iterations++;
274 : }
275 :
276 : /** if we exceed the max iterations, we could still be converged
277 : (considering the acceptable multiplier)
278 : */
279 1 : if (!isConverged(r0_square, r_square, /*acceptable=*/true))
280 0 : _state = State::NOT_CONVERGED;
281 :
282 1 : guess = guess_eigen(0);
283 1 : }
284 :
285 : template <bool is_ad>
286 : template <typename R, typename J>
287 : void
288 1 : NestedSolveTempl<is_ad>::nonlinear(NSRealVectorValue & guess,
289 : R & computeResidual,
290 : J & computeJacobian)
291 : {
292 : using V = NestedSolveTempl<is_ad>::DynamicVector;
293 :
294 : // adaptor functor to utilize the Eigen non-linear solver
295 1 : auto functor = make_realvector_adaptor<true>(computeResidual, computeJacobian);
296 1 : Eigen::HybridNonLinearSolver<decltype(functor), NSReal> solver(functor);
297 1 : V guess_eigen(3);
298 1 : guess_eigen << guess(0), guess(1), guess(2);
299 1 : auto status = solver.solve(guess_eigen);
300 1 : guess(0) = guess_eigen(0);
301 1 : guess(1) = guess_eigen(1);
302 1 : guess(2) = guess_eigen(2);
303 :
304 1 : if (status == Eigen::HybridNonLinearSolverSpace::RelativeErrorTooSmall)
305 1 : _state = State::CONVERGED_REL;
306 : else
307 0 : _state = State::NOT_CONVERGED;
308 1 : }
309 :
310 : template <bool is_ad>
311 : template <typename R, typename J>
312 : void
313 : NestedSolveTempl<is_ad>::nonlinearTight(NestedSolveTempl<is_ad>::DynamicVector & guess,
314 : R & computeResidual,
315 : J & computeJacobian)
316 : {
317 : // adaptor functor to utilize the Eigen non-linear solver
318 : auto functor = make_adaptor<false>(computeResidual, computeJacobian);
319 : Eigen::HybridNonLinearSolver<decltype(functor), NSReal> solver(functor);
320 : auto status = solver.solve(guess);
321 :
322 : if (status == Eigen::HybridNonLinearSolverSpace::RelativeErrorTooSmall)
323 : _state = State::CONVERGED_REL;
324 : else
325 : _state = State::NOT_CONVERGED;
326 : }
327 :
328 : template <bool is_ad>
329 : template <typename R, typename J>
330 : void
331 1 : NestedSolveTempl<is_ad>::nonlinearTight(NSReal & guess, R & computeResidual, J & computeJacobian)
332 : {
333 : using V = NestedSolveTempl<is_ad>::DynamicVector;
334 :
335 : // adaptor functor to utilize the Eigen non-linear solver
336 1 : auto functor = make_real_adaptor<false>(computeResidual, computeJacobian);
337 1 : Eigen::HybridNonLinearSolver<decltype(functor), NSReal> solver(functor);
338 1 : V guess_eigen(1);
339 1 : guess_eigen << guess;
340 1 : auto status = solver.solve(guess_eigen);
341 1 : guess = guess_eigen(0);
342 :
343 1 : if (status == Eigen::HybridNonLinearSolverSpace::RelativeErrorTooSmall)
344 1 : _state = State::CONVERGED_REL;
345 : else
346 0 : _state = State::NOT_CONVERGED;
347 1 : }
348 :
349 : template <bool is_ad>
350 : template <typename R, typename J>
351 : void
352 : NestedSolveTempl<is_ad>::nonlinearTight(NSRealVectorValue & guess,
353 : R & computeResidual,
354 : J & computeJacobian)
355 : {
356 : using V = NestedSolveTempl<is_ad>::DynamicVector;
357 :
358 : // adaptor functor to utilize the Eigen non-linear solver
359 : auto functor = make_realvector_adaptor<false>(computeResidual, computeJacobian);
360 : Eigen::HybridNonLinearSolver<decltype(functor), NSReal> solver(functor);
361 : V guess_eigen(3);
362 : guess_eigen << guess(0), guess(1), guess(2);
363 : auto status = solver.solve(guess_eigen);
364 : guess(0) = guess_eigen(0);
365 : guess(1) = guess_eigen(1);
366 : guess(2) = guess_eigen(2);
367 :
368 : if (status == Eigen::HybridNonLinearSolverSpace::RelativeErrorTooSmall)
369 : _state = State::CONVERGED_REL;
370 : else
371 : _state = State::NOT_CONVERGED;
372 : }
373 :
374 : template <bool is_ad>
375 : template <typename R, typename J, typename B>
376 : void
377 : NestedSolveTempl<is_ad>::nonlinearBounded(NSReal & guess,
378 : R & computeResidual,
379 : J & computeJacobian,
380 : B & computeBounds)
381 : {
382 : auto functor = make_real_adaptor<false>(computeResidual, computeJacobian);
383 : Eigen::HybridNonLinearSolver<decltype(functor), NSReal> solver(functor);
384 :
385 : NestedSolveTempl<is_ad>::DynamicVector guess_eigen(1);
386 : guess_eigen << guess;
387 :
388 : auto status = solver.solveInit(guess_eigen);
389 : if (status == Eigen::HybridNonLinearSolverSpace::ImproperInputParameters)
390 : return;
391 :
392 : bool bounds_hit = false;
393 : while (status == Eigen::HybridNonLinearSolverSpace::Running)
394 : {
395 : status = solver.solveOneStep(guess_eigen);
396 : const std::pair<Real, Real> & bounds = computeBounds();
397 :
398 : // Snap to bounds. Terminate solve if we snap twice consecutively.
399 : if (guess_eigen(0) < bounds.first || guess_eigen(0) > bounds.second)
400 : {
401 : if (guess_eigen(0) < bounds.first)
402 : guess_eigen(0) = bounds.first;
403 : else
404 : guess_eigen(0) = bounds.second;
405 :
406 : if (bounds_hit)
407 : {
408 : _state = State::CONVERGED_BOUNDS;
409 : return;
410 : }
411 :
412 : bounds_hit = true;
413 : }
414 : }
415 :
416 : if (status == Eigen::HybridNonLinearSolverSpace::RelativeErrorTooSmall)
417 : _state = State::CONVERGED_REL;
418 : else
419 : _state = State::NOT_CONVERGED;
420 :
421 : guess = guess_eigen(0);
422 : }
423 :
424 : template <bool is_ad>
425 : template <typename V, typename T>
426 : void
427 7626 : NestedSolveTempl<is_ad>::nonlinear(V & guess, T && compute)
428 : {
429 4 : V delta;
430 4 : V residual;
431 4 : CorrespondingJacobian<V> jacobian;
432 7626 : sizeItems(guess, residual, jacobian);
433 :
434 7626 : _n_iterations = 0;
435 7626 : compute(guess, residual, jacobian);
436 :
437 : // compute first residual norm for relative convergence checks
438 7626 : auto r0_square = normSquare(residual);
439 7626 : if (r0_square == 0)
440 : {
441 0 : _state = State::EXACT_GUESS;
442 7622 : return;
443 : }
444 7626 : auto r_square = r0_square;
445 :
446 : // perform non-linear iterations
447 121301 : while (_n_iterations < _max_iterations)
448 : {
449 : // check convergence
450 121300 : if (_n_iterations >= _min_iterations && isConverged(r0_square, r_square, /*acceptable=*/false))
451 7624 : return;
452 :
453 : // solve and apply next increment
454 113676 : linear(jacobian, delta, residual);
455 :
456 : // Check if step size is smaller than the floating point tolerance
457 113676 : if (isRelSmall(delta, guess, _delta_thresh))
458 : {
459 1 : _state = State::CONVERGED_XTOL;
460 1 : return;
461 : }
462 :
463 113675 : guess -= delta;
464 113675 : _n_iterations++;
465 :
466 : // compute residual and jacobian for the next iteration
467 113675 : compute(guess, residual, jacobian);
468 :
469 113675 : r_square = normSquare(residual);
470 : }
471 :
472 : /** if we exceed the max iterations, we could still be converged
473 : (considering the acceptable multiplier)
474 : */
475 1 : if (!isConverged(r0_square, r_square, /*acceptable=*/true))
476 1 : _state = State::NOT_CONVERGED;
477 7 : }
478 :
479 : template <bool is_ad>
480 : template <typename V, typename T, typename C>
481 : void
482 : NestedSolveTempl<is_ad>::nonlinearDamped(V & guess, T && compute, C && computeCondition)
483 : {
484 : V delta;
485 : V residual;
486 : CorrespondingJacobian<V> jacobian;
487 : sizeItems(guess, residual, jacobian);
488 :
489 : _n_iterations = 0;
490 : compute(guess, residual, jacobian);
491 :
492 : // compute first residual norm for relative convergence checks
493 : auto r0_square = normSquare(residual);
494 : if (r0_square == 0)
495 : {
496 : _state = State::EXACT_GUESS;
497 : return;
498 : }
499 : auto r_square = r0_square;
500 :
501 : // perform non-linear iterations
502 : while (_n_iterations < _max_iterations)
503 : {
504 : // check convergence
505 : if (_n_iterations >= _min_iterations && isConverged(r0_square, r_square, /*acceptable=*/false))
506 : return;
507 :
508 : // solve and apply next increment
509 : linear(jacobian, delta, residual);
510 :
511 : // Check if step size is smaller than the floating point tolerance
512 : if (isRelSmall(delta, guess, _delta_thresh))
513 : {
514 : _state = State::CONVERGED_XTOL;
515 : return;
516 : }
517 :
518 : Real alpha = 1;
519 : unsigned int damping_iterations = 0;
520 : auto prev_guess = guess;
521 : guess -= delta;
522 :
523 : while (!computeCondition(guess) && (damping_iterations < _max_damping_iterations))
524 : {
525 : alpha *= _damping_factor;
526 : guess = prev_guess - alpha * delta;
527 : damping_iterations += 1;
528 : }
529 : _n_iterations++;
530 :
531 : // compute residual and jacobian for the next iteration
532 : compute(guess, residual, jacobian);
533 :
534 : r_square = normSquare(residual);
535 : }
536 :
537 : // if we exceed the max iterations, we could still be converged (considering the acceptable
538 : // multiplier)
539 : if (!isConverged(r0_square, r_square, /*acceptable=*/true))
540 : _state = State::NOT_CONVERGED;
541 : }
542 :
543 : template <bool is_ad>
544 : template <typename J, typename V>
545 : void
546 48 : NestedSolveTempl<is_ad>::linear(const J & A, V & x, const V & b) const
547 : {
548 : // we could make the linear solve method configurable here
549 48 : x = A.partialPivLu().solve(b);
550 48 : }
551 :
552 : template <bool is_ad>
553 : template <typename V>
554 : Real
555 90 : NestedSolveTempl<is_ad>::normSquare(const V & v)
556 : {
557 : // we currently cannot get the raw_value of an AD Eigen Matrix, so we loop over components
558 90 : Real norm_square = 0.0;
559 280 : for (const auto i : index_range(v))
560 190 : norm_square += Utility::pow<2>(MetaPhysicL::raw_value(v.data()[i]));
561 90 : return norm_square;
562 : }
563 :
564 : /**
565 : * Internal use namespace that contains template helpers to map each residual
566 : * type to its corresponding Jacobian type.
567 : */
568 : namespace NestedSolveInternal
569 : {
570 :
571 : template <>
572 : struct CorrespondingJacobianTempl<Real>
573 : {
574 : using type = Real;
575 : };
576 :
577 : template <>
578 : struct CorrespondingJacobianTempl<ADReal>
579 : {
580 : using type = ADReal;
581 : };
582 :
583 : template <>
584 : struct CorrespondingJacobianTempl<RealVectorValue>
585 : {
586 : using type = RankTwoTensor;
587 : };
588 :
589 : template <>
590 : struct CorrespondingJacobianTempl<ADRealVectorValue>
591 : {
592 : using type = ADRankTwoTensor;
593 : };
594 :
595 : template <>
596 : struct CorrespondingJacobianTempl<typename NestedSolveTempl<false>::DynamicVector>
597 : {
598 : using type = NestedSolveTempl<false>::DynamicMatrix;
599 : };
600 :
601 : template <>
602 : struct CorrespondingJacobianTempl<typename NestedSolveTempl<true>::DynamicVector>
603 : {
604 : using type = NestedSolveTempl<true>::DynamicMatrix;
605 : };
606 :
607 : template <int N>
608 : struct CorrespondingJacobianTempl<typename Eigen::Matrix<Real, N, 1>>
609 : {
610 : using type = Eigen::Matrix<Real, N, N>;
611 : };
612 :
613 : template <int N>
614 : struct CorrespondingJacobianTempl<typename Eigen::Matrix<ADReal, N, 1>>
615 : {
616 : using type = Eigen::Matrix<ADReal, N, N>;
617 : };
618 :
619 : /**
620 : * Adaptor functor for Eigen::Matrix based residual and Jacobian types. No type conversion
621 : * required.
622 : */
623 : template <bool is_ad, typename R, typename J, bool store_residual_norm>
624 : class DynamicMatrixEigenAdaptorFunctor
625 : {
626 : using V = typename NestedSolveTempl<is_ad>::DynamicVector;
627 :
628 : public:
629 1 : DynamicMatrixEigenAdaptorFunctor(R & residual_lambda, J & jacobian_lambda)
630 1 : : _residual_lambda(residual_lambda), _jacobian_lambda(jacobian_lambda)
631 : {
632 1 : }
633 29 : int operator()(V & guess, V & residual)
634 : {
635 29 : _residual_lambda(guess, residual);
636 : if constexpr (store_residual_norm)
637 29 : _residual_norm = NestedSolveTempl<is_ad>::normSquare(residual);
638 29 : return 0;
639 : }
640 1 : int df(V & guess, typename NestedSolveTempl<is_ad>::template CorrespondingJacobian<V> & jacobian)
641 : {
642 1 : _jacobian_lambda(guess, jacobian);
643 1 : return 0;
644 : }
645 :
646 : const Real & getResidualNorm()
647 : {
648 : if constexpr (store_residual_norm)
649 : return _residual_norm;
650 : }
651 :
652 : private:
653 : R & _residual_lambda;
654 : J & _jacobian_lambda;
655 :
656 : Real _residual_norm;
657 : };
658 :
659 : /**
660 : * Adaptor functor for scalar Real valued residual and Jacobian types. Picks the single scalar
661 : * component.
662 : */
663 : template <bool is_ad, typename R, typename J, bool store_residual_norm>
664 : class RealEigenAdaptorFunctor
665 : {
666 : using V = typename NestedSolveTempl<is_ad>::DynamicVector;
667 :
668 : public:
669 2 : RealEigenAdaptorFunctor(R & residual_lambda, J & jacobian_lambda)
670 2 : : _residual_lambda(residual_lambda), _jacobian_lambda(jacobian_lambda)
671 : {
672 2 : }
673 22 : int operator()(V & guess, V & residual)
674 : {
675 22 : _residual_lambda(guess(0), residual(0));
676 : if constexpr (store_residual_norm)
677 11 : _residual_norm = NestedSolveTempl<is_ad>::normSquare(residual(0));
678 22 : return 0;
679 : }
680 2 : int df(V & guess, typename NestedSolveTempl<is_ad>::template CorrespondingJacobian<V> & jacobian)
681 : {
682 2 : _jacobian_lambda(guess(0), jacobian(0, 0));
683 2 : return 0;
684 : }
685 :
686 1 : const Real & getResidualNorm()
687 : {
688 : if constexpr (store_residual_norm)
689 1 : return _residual_norm;
690 : }
691 :
692 : private:
693 : R & _residual_lambda;
694 : J & _jacobian_lambda;
695 :
696 : Real _residual_norm;
697 : };
698 :
699 : /**
700 : * Adaptor functor for MOOSE RealVectorValue/RankTwoTensor residual and Jacobian types. Performs
701 : * type conversion.
702 : */
703 : template <bool is_ad, typename R, typename J, bool store_residual_norm>
704 : class RealVectorEigenAdaptorFunctor
705 : {
706 : using V = typename NestedSolveTempl<is_ad>::DynamicVector;
707 :
708 : public:
709 1 : RealVectorEigenAdaptorFunctor(R & residual_lambda, J & jacobian_lambda)
710 1 : : _residual_lambda(residual_lambda), _jacobian_lambda(jacobian_lambda)
711 : {
712 1 : }
713 10 : int operator()(V & guess, V & residual)
714 : {
715 10 : typename NestedSolveTempl<is_ad>::NSRealVectorValue guess_moose(guess(0), guess(1), guess(2));
716 10 : typename NestedSolveTempl<is_ad>::NSRealVectorValue residual_moose;
717 10 : _residual_lambda(guess_moose, residual_moose);
718 10 : residual(0) = residual_moose(0);
719 10 : residual(1) = residual_moose(1);
720 10 : residual(2) = residual_moose(2);
721 : if constexpr (store_residual_norm)
722 10 : _residual_norm = NestedSolveTempl<is_ad>::normSquare(residual);
723 10 : return 0;
724 : }
725 1 : int df(V & guess, typename NestedSolveTempl<is_ad>::template CorrespondingJacobian<V> & jacobian)
726 : {
727 1 : typename NestedSolveTempl<is_ad>::NSRealVectorValue guess_moose(guess(0), guess(1), guess(2));
728 1 : typename NestedSolveTempl<is_ad>::NSRankTwoTensor jacobian_moose;
729 1 : _jacobian_lambda(guess_moose, jacobian_moose);
730 1 : jacobian =
731 1 : Eigen::Map<typename NestedSolveTempl<is_ad>::DynamicMatrix>(&jacobian_moose(0, 0), 3, 3);
732 1 : return 0;
733 1 : }
734 :
735 : const Real & getResidualNorm()
736 : {
737 : if constexpr (store_residual_norm)
738 : return _residual_norm;
739 : }
740 :
741 : private:
742 : R & _residual_lambda;
743 : J & _jacobian_lambda;
744 :
745 : Real _residual_norm;
746 : };
747 :
748 : } // namespace NestedSolveInternal
749 :
750 : template <bool is_ad>
751 : template <bool store_residual_norm, typename R, typename J>
752 : auto
753 1 : NestedSolveTempl<is_ad>::make_adaptor(R & residual_lambda, J & jacobian_lambda)
754 : {
755 : return NestedSolveInternal::DynamicMatrixEigenAdaptorFunctor<is_ad, R, J, store_residual_norm>(
756 1 : residual_lambda, jacobian_lambda);
757 : }
758 : template <bool is_ad>
759 : template <bool store_residual_norm, typename R, typename J>
760 : auto
761 2 : NestedSolveTempl<is_ad>::make_real_adaptor(R & residual_lambda, J & jacobian_lambda)
762 : {
763 : return NestedSolveInternal::RealEigenAdaptorFunctor<is_ad, R, J, store_residual_norm>(
764 2 : residual_lambda, jacobian_lambda);
765 : }
766 :
767 : template <bool is_ad>
768 : template <bool store_residual_norm, typename R, typename J>
769 : auto
770 1 : NestedSolveTempl<is_ad>::make_realvector_adaptor(R & residual_lambda, J & jacobian_lambda)
771 : {
772 : return NestedSolveInternal::RealVectorEigenAdaptorFunctor<is_ad, R, J, store_residual_norm>(
773 1 : residual_lambda, jacobian_lambda);
774 : }
775 :
776 : // lambda to check for convergence and set the state accordingly
777 : template <bool is_ad>
778 : bool
779 98425 : NestedSolveTempl<is_ad>::isConverged(const Real r0_square, Real r_square, bool acceptable)
780 : {
781 98425 : if (acceptable)
782 2 : r_square /= _acceptable_multiplier * _acceptable_multiplier;
783 98425 : if (r_square < _absolute_tolerance_square)
784 : {
785 7624 : _state = acceptable ? State::CONVERGED_ACCEPTABLE_ABS : State::CONVERGED_ABS;
786 7624 : return true;
787 : }
788 90801 : if (r_square / r0_square < _relative_tolerance_square)
789 : {
790 1 : _state = acceptable ? State::CONVERGED_ACCEPTABLE_REL : State::CONVERGED_REL;
791 1 : return true;
792 : }
793 90800 : return false;
794 : }
795 :
796 : typedef NestedSolveTempl<false> NestedSolve;
797 : typedef NestedSolveTempl<true> ADNestedSolve;
|
