Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/diff_system.h"
20 : #include "libmesh/dof_map.h"
21 : #include "libmesh/libmesh_logging.h"
22 : #include "libmesh/linear_solver.h"
23 : #include "libmesh/newton_solver.h"
24 : #include "libmesh/numeric_vector.h"
25 : #include "libmesh/sparse_matrix.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 : // SIGN from Numerical Recipes
31 : template <typename T>
32 : inline
33 124 : T SIGN(T a, T b)
34 : {
35 3597 : return b >= 0 ? std::abs(a) : -std::abs(a);
36 : }
37 :
38 13199 : Real NewtonSolver::line_search(Real tol,
39 : Real last_residual,
40 : Real & current_residual,
41 : NumericVector<Number> & newton_iterate,
42 : const NumericVector<Number> & linear_solution)
43 : {
44 : // Take a full step if we got a residual reduction or if we
45 : // aren't substepping
46 13277 : if ((current_residual < last_residual) ||
47 2834 : (!require_residual_reduction &&
48 0 : (!require_finite_residual || !libmesh_isnan(current_residual))))
49 286 : return 1.;
50 :
51 : // The residual vector
52 2834 : NumericVector<Number> & rhs = *(_system.rhs);
53 :
54 78 : Real ax = 0.; // First abscissa, don't take negative steps
55 78 : Real cx = 1.; // Second abscissa, don't extrapolate steps
56 :
57 : // Find bx, a step length that gives lower residual than ax or cx
58 78 : Real bx = 1.;
59 :
60 5824 : while (libmesh_isnan(current_residual) ||
61 2912 : (current_residual > last_residual &&
62 2834 : require_residual_reduction))
63 : {
64 : // Reduce step size to 1/2, 1/4, etc.
65 : Real substepdivision;
66 2834 : if (brent_line_search && !libmesh_isnan(current_residual))
67 : {
68 2834 : substepdivision = std::min(0.5, last_residual/current_residual);
69 2834 : substepdivision = std::max(substepdivision, tol*2.);
70 : }
71 : else
72 0 : substepdivision = 0.5;
73 :
74 2834 : newton_iterate.add (bx * (1.-substepdivision),
75 156 : linear_solution);
76 2834 : newton_iterate.close();
77 2834 : bx *= substepdivision;
78 2834 : if (verbose)
79 78 : libMesh::out << " Shrinking Newton step to "
80 78 : << bx << std::endl;
81 :
82 : // We may need to localize a parallel solution
83 2834 : _system.update();
84 :
85 : // Check residual with fractional Newton step
86 2834 : _system.assembly(true, false, !this->_exact_constraint_enforcement);
87 :
88 2834 : rhs.close();
89 2834 : current_residual = rhs.l2_norm();
90 2834 : if (verbose)
91 78 : libMesh::out << " Current Residual: "
92 156 : << current_residual << std::endl;
93 :
94 2834 : if (bx/2. < minsteplength &&
95 0 : (libmesh_isnan(current_residual) ||
96 0 : (current_residual > last_residual)))
97 : {
98 0 : libMesh::out << "Inexact Newton step FAILED at step "
99 0 : << _outer_iterations << std::endl;
100 :
101 0 : if (!continue_after_backtrack_failure)
102 : {
103 0 : libmesh_convergence_failure();
104 : }
105 : else
106 : {
107 0 : libMesh::out << "Continuing anyway ..." << std::endl;
108 0 : _solve_result = DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
109 0 : return bx;
110 : }
111 : }
112 : } // end while (current_residual > last_residual)
113 :
114 : // Now return that reduced-residual step, or use Brent's method to
115 : // find a more optimal step.
116 :
117 2834 : if (!brent_line_search)
118 0 : return bx;
119 :
120 : // Brent's method adapted from Numerical Recipes in C, ch. 10.2
121 78 : Real e = 0.;
122 :
123 78 : Real x = bx, w = bx, v = bx;
124 :
125 : // Residuals at bx
126 78 : Real fx = current_residual,
127 78 : fw = current_residual,
128 78 : fv = current_residual;
129 :
130 : // Max iterations for Brent's method loop
131 78 : const unsigned int max_i = 20;
132 :
133 : // for golden ratio steps
134 78 : const Real golden_ratio = 1.-(std::sqrt(5.)-1.)/2.;
135 :
136 22841 : for (unsigned int i=1; i <= max_i; i++)
137 : {
138 22841 : Real xm = (ax+cx)*0.5;
139 22841 : Real tol1 = tol * std::abs(x) + tol*tol;
140 22841 : Real tol2 = 2.0 * tol1;
141 :
142 : // Test if we're done
143 23469 : if (std::abs(x-xm) <= (tol2 - 0.5 * (cx - ax)))
144 78 : return x;
145 :
146 : Real d;
147 :
148 : // Construct a parabolic fit
149 20007 : if (std::abs(e) > tol1)
150 : {
151 17173 : Real r = (x-w)*(fx-fv);
152 17173 : Real q = (x-v)*(fx-fw);
153 17173 : Real p = (x-v)*q-(x-w)*r;
154 17173 : q = 2. * (q-r);
155 17173 : if (q > 0.)
156 6471 : p = -p;
157 : else
158 290 : q = std::abs(q);
159 30724 : if (std::abs(p) >= std::abs(0.5*q*e) ||
160 17567 : p <= q * (ax-x) ||
161 14339 : p >= q * (cx-x))
162 : {
163 : // Take a golden section step
164 2834 : e = x >= xm ? ax-x : cx-x;
165 2834 : d = golden_ratio * e;
166 : }
167 : else
168 : {
169 : // Take a parabolic fit step
170 14339 : d = p/q;
171 14339 : if (x+d-ax < tol2 || cx-(x+d) < tol2)
172 2863 : d = SIGN(tol1, xm - x);
173 : }
174 : }
175 : else
176 : {
177 : // Take a golden section step
178 2834 : e = x >= xm ? ax-x : cx-x;
179 2834 : d = golden_ratio * e;
180 : }
181 :
182 21881 : Real u = std::abs(d) >= tol1 ? x+d : x + SIGN(tol1,d);
183 :
184 : // Assemble the residual at the new steplength u
185 20007 : newton_iterate.add (bx - u, linear_solution);
186 20007 : newton_iterate.close();
187 550 : bx = u;
188 20007 : if (verbose)
189 550 : libMesh::out << " Shrinking Newton step to "
190 550 : << bx << std::endl;
191 :
192 : // We may need to localize a parallel solution
193 20007 : _system.update();
194 20007 : _system.assembly(true, false, !this->_exact_constraint_enforcement);
195 :
196 20007 : rhs.close();
197 20007 : Real fu = current_residual = rhs.l2_norm();
198 20007 : if (verbose)
199 550 : libMesh::out << " Current Residual: "
200 550 : << fu << std::endl;
201 :
202 20007 : if (fu <= fx)
203 : {
204 10297 : if (u >= x)
205 122 : ax = x;
206 : else
207 162 : cx = x;
208 284 : v = w; w = x; x = u;
209 284 : fv = fw; fw = fx; fx = fu;
210 : }
211 : else
212 : {
213 9710 : if (u < x)
214 104 : ax = u;
215 : else
216 162 : cx = u;
217 9710 : if (fu <= fw || w == x)
218 : {
219 212 : v = w; w = u;
220 212 : fv = fw; fw = fu;
221 : }
222 1965 : else if (fu <= fv || v == x || v == w)
223 : {
224 54 : v = u;
225 54 : fv = fu;
226 : }
227 : }
228 : }
229 :
230 0 : if (!quiet)
231 0 : libMesh::out << "Warning! Too many iterations used in Brent line search!"
232 0 : << std::endl;
233 0 : return bx;
234 : }
235 :
236 :
237 4975 : NewtonSolver::NewtonSolver (sys_type & s)
238 : : Parent(s),
239 4679 : require_residual_reduction(true),
240 4679 : require_finite_residual(true),
241 4679 : brent_line_search(true),
242 4679 : track_linear_convergence(false),
243 4679 : minsteplength(1e-5),
244 4679 : linear_tolerance_multiplier(1e-3),
245 4975 : _linear_solver(LinearSolver<Number>::build(s.comm()))
246 : {
247 4975 : }
248 :
249 :
250 :
251 9358 : NewtonSolver::~NewtonSolver () = default;
252 :
253 :
254 :
255 4975 : void NewtonSolver::init()
256 : {
257 4975 : Parent::init();
258 :
259 9802 : if (libMesh::on_command_line("--solver-system-names"))
260 0 : _linear_solver->init((_system.name()+"_").c_str());
261 : else
262 4975 : _linear_solver->init();
263 :
264 4975 : _linear_solver->init_names(_system);
265 4975 : }
266 :
267 :
268 :
269 5236 : void NewtonSolver::reinit()
270 : {
271 5236 : Parent::reinit();
272 :
273 5236 : _linear_solver->clear();
274 :
275 5236 : _linear_solver->init_names(_system);
276 5236 : }
277 :
278 :
279 :
280 29803 : unsigned int NewtonSolver::solve()
281 : {
282 1764 : LOG_SCOPE("solve()", "NewtonSolver");
283 :
284 : // Reset any prior solve result
285 29803 : _solve_result = INVALID_SOLVE_RESULT;
286 :
287 29803 : NumericVector<Number> & newton_iterate = *(_system.solution);
288 :
289 29803 : std::unique_ptr<NumericVector<Number>> linear_solution_ptr = newton_iterate.zero_clone();
290 882 : NumericVector<Number> & linear_solution = *linear_solution_ptr;
291 29803 : NumericVector<Number> & rhs = *(_system.rhs);
292 :
293 29803 : newton_iterate.close();
294 29803 : linear_solution.close();
295 29803 : rhs.close();
296 :
297 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
298 29803 : if (this->_exact_constraint_enforcement)
299 27805 : _system.get_dof_map().enforce_constraints_exactly(_system);
300 : #endif
301 :
302 29803 : SparseMatrix<Number> & matrix = *(_system.matrix);
303 :
304 : // Set starting linear tolerance
305 29803 : double current_linear_tolerance = initial_linear_tolerance;
306 :
307 : // Start counting our linear solver steps
308 29803 : _inner_iterations = 0;
309 :
310 : // Now we begin the nonlinear loop
311 41937 : for (_outer_iterations=0; _outer_iterations<max_nonlinear_iterations;
312 12134 : ++_outer_iterations)
313 : {
314 : // We may need to localize a parallel solution
315 41937 : _system.update();
316 :
317 41937 : if (verbose)
318 338 : libMesh::out << "Assembling the System" << std::endl;
319 :
320 41937 : _system.assembly(true, true, !this->_exact_constraint_enforcement);
321 41937 : rhs.close();
322 41937 : Real current_residual = rhs.l2_norm();
323 :
324 40721 : if (libmesh_isnan(current_residual))
325 : {
326 0 : libMesh::out << " Nonlinear solver DIVERGED at step "
327 0 : << _outer_iterations
328 0 : << " with residual Not-a-Number"
329 0 : << std::endl;
330 0 : libmesh_convergence_failure();
331 0 : continue;
332 : }
333 :
334 41937 : if (current_residual <= absolute_residual_tolerance)
335 : {
336 0 : if (verbose)
337 0 : libMesh::out << "Linear solve unnecessary; residual "
338 0 : << current_residual
339 0 : << " meets absolute tolerance "
340 0 : << absolute_residual_tolerance
341 0 : << std::endl;
342 :
343 : // We're not doing a solve, but other code may reuse this
344 : // matrix.
345 0 : matrix.close();
346 :
347 0 : _solve_result |= CONVERGED_ABSOLUTE_RESIDUAL;
348 0 : if (current_residual == 0)
349 : {
350 0 : if (relative_residual_tolerance > 0)
351 0 : _solve_result |= CONVERGED_RELATIVE_RESIDUAL;
352 0 : if (absolute_step_tolerance > 0)
353 0 : _solve_result |= CONVERGED_ABSOLUTE_STEP;
354 0 : if (relative_step_tolerance > 0)
355 0 : _solve_result |= CONVERGED_RELATIVE_STEP;
356 : }
357 :
358 29803 : break;
359 : }
360 :
361 : // Prepare to take incomplete steps
362 1216 : Real last_residual = current_residual;
363 :
364 41937 : max_residual_norm = std::max (current_residual,
365 41937 : max_residual_norm);
366 :
367 : // Compute the l2 norm of the whole solution
368 41937 : Real norm_total = newton_iterate.l2_norm();
369 41937 : max_solution_norm = std::max(max_solution_norm, norm_total);
370 :
371 41937 : if (verbose)
372 338 : libMesh::out << "Nonlinear Residual: "
373 676 : << current_residual << std::endl;
374 :
375 : // Make sure our linear tolerance is low enough
376 41937 : current_linear_tolerance =
377 41937 : double(std::min (current_linear_tolerance,
378 41937 : current_residual * linear_tolerance_multiplier));
379 :
380 : // But don't let it be too small
381 41937 : if (current_linear_tolerance < minimum_linear_tolerance)
382 : {
383 16982 : current_linear_tolerance = minimum_linear_tolerance;
384 : }
385 :
386 : // If starting the nonlinear solve with a really good initial guess, we dont want to set an absurd linear tolerance
387 41937 : current_linear_tolerance =
388 41937 : double(std::max(current_linear_tolerance,
389 41937 : absolute_residual_tolerance / current_residual
390 41937 : / 10.0));
391 :
392 : // At this point newton_iterate is the current guess, and
393 : // linear_solution is now about to become the NEGATIVE of the next
394 : // Newton step.
395 :
396 : // Our best initial guess for the linear_solution is zero!
397 41937 : linear_solution.zero();
398 :
399 41937 : if (verbose)
400 338 : libMesh::out << "Linear solve starting, tolerance "
401 338 : << current_linear_tolerance << std::endl;
402 :
403 : // Solve the linear system.
404 : const std::pair<unsigned int, Real> rval =
405 43153 : _linear_solver->solve (matrix, _system.request_matrix("Preconditioner"),
406 : linear_solution, rhs, current_linear_tolerance,
407 2432 : max_linear_iterations);
408 :
409 41937 : if (track_linear_convergence)
410 : {
411 0 : LinearConvergenceReason linear_c_reason = _linear_solver->get_converged_reason();
412 :
413 : // Check if something went wrong during the linear solve
414 0 : if (linear_c_reason < 0)
415 : {
416 : // The linear solver failed somehow
417 0 : _solve_result |= DiffSolver::DIVERGED_LINEAR_SOLVER_FAILURE;
418 : // Print a message
419 0 : libMesh::out << "Linear solver failed during Newton step, dropping out."
420 0 : << std::endl;
421 0 : break;
422 : }
423 : }
424 :
425 : // We may need to localize a parallel solution
426 41937 : _system.update ();
427 : // The linear solver may not have fit our constraints exactly
428 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
429 41937 : if (this->_exact_constraint_enforcement)
430 39137 : _system.get_dof_map().enforce_constraints_exactly
431 39137 : (_system, &linear_solution, /* homogeneous = */ true);
432 : #endif
433 :
434 41937 : const unsigned int linear_steps = rval.first;
435 1216 : libmesh_assert_less_equal (linear_steps, max_linear_iterations);
436 41937 : _inner_iterations += linear_steps;
437 :
438 41937 : const bool linear_solve_finished =
439 41937 : !(linear_steps == max_linear_iterations);
440 :
441 41937 : if (verbose)
442 338 : libMesh::out << "Linear solve finished, step " << linear_steps
443 676 : << ", residual " << rval.second
444 338 : << std::endl;
445 :
446 : // Compute the l2 norm of the nonlinear update
447 41937 : Real norm_delta = linear_solution.l2_norm();
448 :
449 41937 : if (verbose)
450 338 : libMesh::out << "Trying full Newton step" << std::endl;
451 : // Take a full Newton step
452 41937 : newton_iterate.add (-1., linear_solution);
453 41937 : newton_iterate.close();
454 :
455 41937 : if (this->linear_solution_monitor.get())
456 : {
457 : // Compute the l2 norm of the whole solution
458 0 : norm_total = newton_iterate.l2_norm();
459 0 : rhs.close();
460 0 : (*this->linear_solution_monitor)(linear_solution, norm_delta,
461 : newton_iterate, norm_total,
462 0 : rhs, rhs.l2_norm(), _outer_iterations);
463 : }
464 :
465 : // Check residual with full Newton step, if that's useful for determining
466 : // whether to line search, whether to quit early, or whether to die after
467 : // hitting our max iteration count
468 41937 : if (this->require_residual_reduction ||
469 650 : this->require_finite_residual ||
470 0 : _outer_iterations+1 < max_nonlinear_iterations ||
471 0 : !continue_after_max_iterations)
472 : {
473 41937 : _system.update ();
474 41937 : _system.assembly(true, false, !this->_exact_constraint_enforcement);
475 :
476 41937 : rhs.close();
477 41937 : current_residual = rhs.l2_norm();
478 41937 : if (verbose)
479 338 : libMesh::out << " Current Residual: "
480 676 : << current_residual << std::endl;
481 :
482 : // don't fiddle around if we've already converged
483 41937 : if (test_convergence(current_residual, norm_delta,
484 43147 : linear_solve_finished &&
485 41913 : current_residual <= last_residual))
486 : {
487 : // Make sure we have an up to date max_solution_norm if
488 : // we're going to be verbose about it here. We don't
489 : // want to update it earlier because we don't want an
490 : // excessive full Newton step to explode it if we're
491 : // just going to line search that step back.
492 28738 : if (!quiet && verbose)
493 : {
494 5632 : norm_total = newton_iterate.l2_norm();
495 5632 : max_solution_norm = std::max(max_solution_norm, norm_total);
496 : }
497 :
498 28738 : if (!quiet)
499 19082 : print_convergence(_outer_iterations, current_residual,
500 19662 : norm_delta, linear_solve_finished &&
501 19082 : current_residual <= last_residual);
502 28738 : _outer_iterations++;
503 28738 : break; // out of _outer_iterations for loop
504 : }
505 : }
506 :
507 : // since we're not converged, backtrack if necessary
508 : Real steplength =
509 13199 : this->line_search(std::sqrt(TOLERANCE),
510 : last_residual, current_residual,
511 : newton_iterate, linear_solution);
512 13199 : norm_delta *= steplength;
513 :
514 : // Check to see if backtracking failed,
515 : // and break out of the nonlinear loop if so...
516 13199 : if (_solve_result == DiffSolver::DIVERGED_BACKTRACKING_FAILURE)
517 : {
518 0 : _outer_iterations++;
519 0 : break; // out of _outer_iterations for loop
520 : }
521 :
522 13199 : if (_outer_iterations + 1 >= max_nonlinear_iterations)
523 : {
524 0 : libMesh::out << " Nonlinear solver reached maximum step "
525 0 : << max_nonlinear_iterations << ", latest evaluated residual "
526 0 : << current_residual << std::endl;
527 0 : if (continue_after_max_iterations)
528 : {
529 0 : _solve_result = DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
530 0 : libMesh::out << " Continuing..." << std::endl;
531 : }
532 : else
533 : {
534 0 : libmesh_convergence_failure();
535 : }
536 0 : continue;
537 : }
538 :
539 : // Compute the l2 norm of the whole solution
540 13199 : norm_total = newton_iterate.l2_norm();
541 :
542 13199 : max_solution_norm = std::max(max_solution_norm, norm_total);
543 :
544 : // Print out information for the
545 : // nonlinear iterations.
546 13199 : if (verbose)
547 142 : libMesh::out << " Nonlinear step: |du|/|u| = "
548 5460 : << norm_delta / norm_total
549 284 : << ", |du| = " << norm_delta
550 142 : << std::endl;
551 :
552 : // Terminate the solution iteration if the difference between
553 : // this iteration and the last is sufficiently small.
554 13199 : if (test_convergence(current_residual, norm_delta / steplength,
555 : linear_solve_finished))
556 : {
557 1065 : if (!quiet)
558 1065 : print_convergence(_outer_iterations, current_residual,
559 : norm_delta / steplength,
560 : linear_solve_finished);
561 1065 : _outer_iterations++;
562 1065 : break; // out of _outer_iterations for loop
563 : }
564 : } // end nonlinear loop
565 :
566 : // The linear solver may not have fit our constraints exactly
567 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
568 29803 : if (this->_exact_constraint_enforcement)
569 27805 : _system.get_dof_map().enforce_constraints_exactly(_system);
570 : #endif
571 :
572 : // We may need to localize a parallel solution
573 29803 : _system.update ();
574 :
575 : // Make sure we are returning something sensible as the
576 : // _solve_result, except in the edge case where we weren't really asked to
577 : // solve.
578 882 : libmesh_assert (_solve_result != DiffSolver::INVALID_SOLVE_RESULT ||
579 : !max_nonlinear_iterations);
580 :
581 31567 : return _solve_result;
582 28039 : }
583 :
584 :
585 :
586 55136 : bool NewtonSolver::test_convergence(Real current_residual,
587 : Real step_norm,
588 : bool linear_solve_finished)
589 : {
590 : // We haven't converged unless we pass a convergence test
591 1580 : bool has_converged = false;
592 :
593 : // Is our absolute residual low enough?
594 55136 : if (current_residual < absolute_residual_tolerance)
595 : {
596 6955 : _solve_result |= CONVERGED_ABSOLUTE_RESIDUAL;
597 196 : has_converged = true;
598 : }
599 :
600 : // Is our relative residual low enough?
601 58296 : if ((current_residual / max_residual_norm) <
602 55136 : relative_residual_tolerance)
603 : {
604 29454 : _solve_result |= CONVERGED_RELATIVE_RESIDUAL;
605 880 : has_converged = true;
606 : }
607 :
608 : // For incomplete linear solves, it's not safe to test step sizes
609 55136 : if (!linear_solve_finished)
610 : {
611 84 : return has_converged;
612 : }
613 :
614 : // Is our absolute Newton step size small enough?
615 52278 : if (step_norm < absolute_step_tolerance)
616 : {
617 0 : _solve_result |= CONVERGED_ABSOLUTE_STEP;
618 0 : has_converged = true;
619 : }
620 :
621 : // Is our relative Newton step size small enough?
622 52278 : if (max_solution_norm &&
623 49483 : (step_norm / max_solution_norm <
624 49483 : relative_step_tolerance))
625 : {
626 1334 : _solve_result |= CONVERGED_RELATIVE_STEP;
627 12 : has_converged = true;
628 : }
629 :
630 1496 : return has_converged;
631 : }
632 :
633 :
634 20147 : void NewtonSolver::print_convergence(unsigned int step_num,
635 : Real current_residual,
636 : Real step_norm,
637 : bool linear_solve_finished)
638 : {
639 : // Is our absolute residual low enough?
640 20147 : if (current_residual < absolute_residual_tolerance)
641 : {
642 12 : libMesh::out << " Nonlinear solver converged, step " << step_num
643 12 : << ", residual " << current_residual
644 12 : << std::endl;
645 : }
646 19730 : else if (absolute_residual_tolerance)
647 : {
648 328 : if (verbose)
649 0 : libMesh::out << " Nonlinear solver current_residual "
650 0 : << current_residual << " > "
651 0 : << (absolute_residual_tolerance) << std::endl;
652 : }
653 :
654 : // Is our relative residual low enough?
655 21367 : if ((current_residual / max_residual_norm) <
656 20147 : relative_residual_tolerance)
657 : {
658 608 : libMesh::out << " Nonlinear solver converged, step " << step_num
659 608 : << ", residual reduction "
660 20406 : << current_residual / max_residual_norm
661 1216 : << " < " << relative_residual_tolerance
662 608 : << std::endl;
663 : }
664 349 : else if (relative_residual_tolerance)
665 : {
666 349 : if (verbose)
667 2 : libMesh::out << " Nonlinear solver relative residual "
668 351 : << (current_residual / max_residual_norm)
669 4 : << " > " << relative_residual_tolerance
670 2 : << std::endl;
671 : }
672 :
673 : // For incomplete linear solves, it's not safe to test step sizes
674 20147 : if (!linear_solve_finished)
675 2 : return;
676 :
677 : // Is our absolute Newton step size small enough?
678 20131 : if (step_norm < absolute_step_tolerance)
679 : {
680 0 : libMesh::out << " Nonlinear solver converged, step " << step_num
681 0 : << ", absolute step size "
682 0 : << step_norm
683 0 : << " < " << absolute_step_tolerance
684 0 : << std::endl;
685 : }
686 20131 : else if (absolute_step_tolerance)
687 : {
688 0 : if (verbose)
689 0 : libMesh::out << " Nonlinear solver absolute step size "
690 0 : << step_norm
691 0 : << " > " << absolute_step_tolerance
692 0 : << std::endl;
693 : }
694 :
695 : // Is our relative Newton step size small enough?
696 20131 : if (max_solution_norm &&
697 19918 : (step_norm / max_solution_norm <
698 19918 : relative_step_tolerance))
699 : {
700 12 : libMesh::out << " Nonlinear solver converged, step " << step_num
701 12 : << ", relative step size "
702 1334 : << (step_norm / max_solution_norm)
703 24 : << " < " << relative_step_tolerance
704 12 : << std::endl;
705 : }
706 18797 : else if (relative_step_tolerance)
707 : {
708 18797 : if (verbose)
709 : {
710 5347 : if (max_solution_norm)
711 212 : libMesh::out << " Nonlinear solver relative step size "
712 5559 : << (step_norm / max_solution_norm)
713 424 : << " > " << relative_step_tolerance
714 212 : << std::endl;
715 : }
716 : }
717 : }
718 :
719 : } // namespace libMesh
|
