newton_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
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 T SIGN(T a, T b)
34 {
35  return b >= 0 ? std::abs(a) : -std::abs(a);
36 }
37 
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  if ((current_residual < last_residual) ||
48  (!require_finite_residual || !libmesh_isnan(current_residual))))
49  return 1.;
50 
51  // The residual vector
53 
54  Real ax = 0.; // First abscissa, don't take negative steps
55  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  Real bx = 1.;
59 
60  while (libmesh_isnan(current_residual) ||
61  (current_residual > last_residual &&
63  {
64  // Reduce step size to 1/2, 1/4, etc.
65  Real substepdivision;
66  if (brent_line_search && !libmesh_isnan(current_residual))
67  {
68  substepdivision = std::min(0.5, last_residual/current_residual);
69  substepdivision = std::max(substepdivision, tol*2.);
70  }
71  else
72  substepdivision = 0.5;
73 
74  newton_iterate.add (bx * (1.-substepdivision),
75  linear_solution);
76  newton_iterate.close();
77  bx *= substepdivision;
78  if (verbose)
79  libMesh::out << " Shrinking Newton step to "
80  << bx << std::endl;
81 
82  // Check residual with fractional Newton step
83  _system.assembly (true, false);
84 
85  rhs.close();
86  current_residual = rhs.l2_norm();
87  if (verbose)
88  libMesh::out << " Current Residual: "
89  << current_residual << std::endl;
90 
91  if (bx/2. < minsteplength &&
92  (libmesh_isnan(current_residual) ||
93  (current_residual > last_residual)))
94  {
95  libMesh::out << "Inexact Newton step FAILED at step "
96  << _outer_iterations << std::endl;
97 
99  {
100  libmesh_convergence_failure();
101  }
102  else
103  {
104  libMesh::out << "Continuing anyway ..." << std::endl;
106  return bx;
107  }
108  }
109  } // end while (current_residual > last_residual)
110 
111  // Now return that reduced-residual step, or use Brent's method to
112  // find a more optimal step.
113 
114  if (!brent_line_search)
115  return bx;
116 
117  // Brent's method adapted from Numerical Recipes in C, ch. 10.2
118  Real e = 0.;
119 
120  Real x = bx, w = bx, v = bx;
121 
122  // Residuals at bx
123  Real fx = current_residual,
124  fw = current_residual,
125  fv = current_residual;
126 
127  // Max iterations for Brent's method loop
128  const unsigned int max_i = 20;
129 
130  // for golden ratio steps
131  const Real golden_ratio = 1.-(std::sqrt(5.)-1.)/2.;
132 
133  for (unsigned int i=1; i <= max_i; i++)
134  {
135  Real xm = (ax+cx)*0.5;
136  Real tol1 = tol * std::abs(x) + tol*tol;
137  Real tol2 = 2.0 * tol1;
138 
139  // Test if we're done
140  if (std::abs(x-xm) <= (tol2 - 0.5 * (cx - ax)))
141  return x;
142 
143  Real d;
144 
145  // Construct a parabolic fit
146  if (std::abs(e) > tol1)
147  {
148  Real r = (x-w)*(fx-fv);
149  Real q = (x-v)*(fx-fw);
150  Real p = (x-v)*q-(x-w)*r;
151  q = 2. * (q-r);
152  if (q > 0.)
153  p = -p;
154  else
155  q = std::abs(q);
156  if (std::abs(p) >= std::abs(0.5*q*e) ||
157  p <= q * (ax-x) ||
158  p >= q * (cx-x))
159  {
160  // Take a golden section step
161  e = x >= xm ? ax-x : cx-x;
162  d = golden_ratio * e;
163  }
164  else
165  {
166  // Take a parabolic fit step
167  d = p/q;
168  if (x+d-ax < tol2 || cx-(x+d) < tol2)
169  d = SIGN(tol1, xm - x);
170  }
171  }
172  else
173  {
174  // Take a golden section step
175  e = x >= xm ? ax-x : cx-x;
176  d = golden_ratio * e;
177  }
178 
179  Real u = std::abs(d) >= tol1 ? x+d : x + SIGN(tol1,d);
180 
181  // Assemble the residual at the new steplength u
182  newton_iterate.add (bx - u, linear_solution);
183  newton_iterate.close();
184  bx = u;
185  if (verbose)
186  libMesh::out << " Shrinking Newton step to "
187  << bx << std::endl;
188 
189  _system.assembly (true, false);
190 
191  rhs.close();
192  Real fu = current_residual = rhs.l2_norm();
193  if (verbose)
194  libMesh::out << " Current Residual: "
195  << fu << std::endl;
196 
197  if (fu <= fx)
198  {
199  if (u >= x)
200  ax = x;
201  else
202  cx = x;
203  v = w; w = x; x = u;
204  fv = fw; fw = fx; fx = fu;
205  }
206  else
207  {
208  if (u < x)
209  ax = u;
210  else
211  cx = u;
212  if (fu <= fw || w == x)
213  {
214  v = w; w = u;
215  fv = fw; fw = fu;
216  }
217  else if (fu <= fv || v == x || v == w)
218  {
219  v = u;
220  fv = fu;
221  }
222  }
223  }
224 
225  if (!quiet)
226  libMesh::out << "Warning! Too many iterations used in Brent line search!"
227  << std::endl;
228  return bx;
229 }
230 
231 
233  : Parent(s),
234  require_residual_reduction(true),
235  require_finite_residual(true),
236  brent_line_search(true),
237  minsteplength(1e-5),
238  linear_tolerance_multiplier(1e-3),
239  linear_solver(LinearSolver<Number>::build(s.comm()))
240 {
241 }
242 
243 
244 
246 {
247 }
248 
249 
250 
252 {
253  Parent::reinit();
254 
255  linear_solver->clear();
256 }
257 
258 
259 
260 unsigned int NewtonSolver::solve()
261 {
262  START_LOG("solve()", "NewtonSolver");
263 
264  // Reset any prior solve result
266 
267  NumericVector<Number> &newton_iterate = *(_system.solution);
268 
269  AutoPtr<NumericVector<Number> > linear_solution_ptr = newton_iterate.zero_clone();
270  NumericVector<Number> &linear_solution = *linear_solution_ptr;
272 
273  newton_iterate.close();
274  linear_solution.close();
275  rhs.close();
276 
277 #ifdef LIBMESH_ENABLE_CONSTRAINTS
279 #endif
280 
281  SparseMatrix<Number> &matrix = *(_system.matrix);
282 
283  // Prepare to take incomplete steps
284  Real last_residual=0.;
285 
286  // Set starting linear tolerance
287  Real current_linear_tolerance = initial_linear_tolerance;
288 
289  // Start counting our linear solver steps
290  _inner_iterations = 0;
291 
292  // Now we begin the nonlinear loop
295  {
296  if (verbose)
297  libMesh::out << "Assembling the System" << std::endl;
298 
299  _system.assembly(true, true);
300  rhs.close();
301  Real current_residual = rhs.l2_norm();
302  last_residual = current_residual;
303 
304  if (libmesh_isnan(current_residual))
305  {
306  libMesh::out << " Nonlinear solver DIVERGED at step "
308  << " with residual Not-a-Number"
309  << std::endl;
310  libmesh_convergence_failure();
311  continue;
312  }
313 
314  max_residual_norm = std::max (current_residual,
316 
317  // Compute the l2 norm of the whole solution
318  Real norm_total = newton_iterate.l2_norm();
319 
321 
322  if (verbose)
323  libMesh::out << "Nonlinear Residual: "
324  << current_residual << std::endl;
325 
326  // Make sure our linear tolerance is low enough
327  current_linear_tolerance = std::min (current_linear_tolerance,
328  current_residual * linear_tolerance_multiplier);
329 
330  // But don't let it be too small
331  if (current_linear_tolerance < minimum_linear_tolerance)
332  {
333  current_linear_tolerance = minimum_linear_tolerance;
334  }
335 
336  // At this point newton_iterate is the current guess, and
337  // linear_solution is now about to become the NEGATIVE of the next
338  // Newton step.
339 
340  // Our best initial guess for the linear_solution is zero!
341  linear_solution.zero();
342 
343  if (verbose)
344  libMesh::out << "Linear solve starting, tolerance "
345  << current_linear_tolerance << std::endl;
346 
347  // Solve the linear system.
348  const std::pair<unsigned int, Real> rval =
349  linear_solver->solve (matrix, _system.request_matrix("Preconditioner"),
350  linear_solution, rhs, current_linear_tolerance,
352 
353  // We may need to localize a parallel solution
354  _system.update ();
355  // The linear solver may not have fit our constraints exactly
356 #ifdef LIBMESH_ENABLE_CONSTRAINTS
358  (_system, &linear_solution, /* homogeneous = */ true);
359 #endif
360 
361  const unsigned int linear_steps = rval.first;
362  libmesh_assert_less_equal (linear_steps, max_linear_iterations);
363  _inner_iterations += linear_steps;
364 
365  const bool linear_solve_finished =
366  !(linear_steps == max_linear_iterations);
367 
368  if (verbose)
369  libMesh::out << "Linear solve finished, step " << linear_steps
370  << ", residual " << rval.second
371  << std::endl;
372 
373  // Compute the l2 norm of the nonlinear update
374  Real norm_delta = linear_solution.l2_norm();
375 
376  if (verbose)
377  libMesh::out << "Trying full Newton step" << std::endl;
378  // Take a full Newton step
379  newton_iterate.add (-1., linear_solution);
380  newton_iterate.close();
381 
382  if (this->linear_solution_monitor.get()) {
383  // Compute the l2 norm of the whole solution
384  norm_total = newton_iterate.l2_norm();
385  rhs.close();
386  (*this->linear_solution_monitor)(
387  linear_solution, norm_delta,
388  newton_iterate, norm_total,
389  rhs, rhs.l2_norm(), _outer_iterations);
390  }
391 
392  // Check residual with full Newton step, if that's useful for determining
393  // whether to line search, whether to quit early, or whether to die after
394  // hitting our max iteration count
395  if (this->require_residual_reduction ||
396  this->require_finite_residual ||
397  _outer_iterations+1 < max_nonlinear_iterations ||
399  {
400  _system.assembly(true, false);
401 
402  rhs.close();
403  current_residual = rhs.l2_norm();
404  if (verbose)
405  libMesh::out << " Current Residual: "
406  << current_residual << std::endl;
407 
408  // don't fiddle around if we've already converged
409  if (test_convergence(current_residual, norm_delta,
410  linear_solve_finished &&
411  current_residual <= last_residual))
412  {
413  if (!quiet)
414  print_convergence(_outer_iterations, current_residual,
415  norm_delta, linear_solve_finished &&
416  current_residual <= last_residual);
418  break; // out of _outer_iterations for loop
419  }
420  }
421 
422  // since we're not converged, backtrack if necessary
423  Real steplength =
424  this->line_search(std::sqrt(TOLERANCE),
425  last_residual, current_residual,
426  newton_iterate, linear_solution);
427  norm_delta *= steplength;
428 
429  // Check to see if backtracking failed,
430  // and break out of the nonlinear loop if so...
432  {
434  break; // out of _outer_iterations for loop
435  }
436 
437  if (_outer_iterations + 1 >= max_nonlinear_iterations)
438  {
439  libMesh::out << " Nonlinear solver reached maximum step "
440  << max_nonlinear_iterations << ", latest evaluated residual "
441  << current_residual << std::endl;
443  {
445  libMesh::out << " Continuing..." << std::endl;
446  }
447  else
448  {
449  libmesh_convergence_failure();
450  }
451  continue;
452  }
453 
454  // Compute the l2 norm of the whole solution
455  norm_total = newton_iterate.l2_norm();
456 
458 
459  // Print out information for the
460  // nonlinear iterations.
461  if (verbose)
462  libMesh::out << " Nonlinear step: |du|/|u| = "
463  << norm_delta / norm_total
464  << ", |du| = " << norm_delta
465  << std::endl;
466 
467  // Terminate the solution iteration if the difference between
468  // this iteration and the last is sufficiently small.
469  if (test_convergence(current_residual, norm_delta / steplength,
470  linear_solve_finished))
471  {
472  if (!quiet)
473  print_convergence(_outer_iterations, current_residual,
474  norm_delta / steplength,
475  linear_solve_finished);
477  break; // out of _outer_iterations for loop
478  }
479  } // end nonlinear loop
480 
481  // The linear solver may not have fit our constraints exactly
482 #ifdef LIBMESH_ENABLE_CONSTRAINTS
484 #endif
485 
486  // We may need to localize a parallel solution
487  _system.update ();
488 
489  STOP_LOG("solve()", "NewtonSolver");
490 
491  // Make sure we are returning something sensible as the
492  // _solve_result, except in the edge case where we weren't really asked to
493  // solve.
495  !max_nonlinear_iterations);
496 
497  return _solve_result;
498 }
499 
500 
501 
502 bool NewtonSolver::test_convergence(Real current_residual,
503  Real step_norm,
504  bool linear_solve_finished)
505 {
506  // We haven't converged unless we pass a convergence test
507  bool has_converged = false;
508 
509  // Is our absolute residual low enough?
510  if (current_residual < absolute_residual_tolerance)
511  {
513  has_converged = true;
514  }
515 
516  // Is our relative residual low enough?
517  if ((current_residual / max_residual_norm) <
519  {
521  has_converged = true;
522  }
523 
524  // For incomplete linear solves, it's not safe to test step sizes
525  if (!linear_solve_finished)
526  {
527  return has_converged;
528  }
529 
530  // Is our absolute Newton step size small enough?
531  if (step_norm < absolute_step_tolerance)
532  {
534  has_converged = true;
535  }
536 
537  // Is our relative Newton step size small enough?
538  if (step_norm / max_solution_norm <
540  {
542  has_converged = true;
543  }
544 
545  return has_converged;
546 }
547 
548 
549 void NewtonSolver::print_convergence(unsigned int step_num,
550  Real current_residual,
551  Real step_norm,
552  bool linear_solve_finished)
553 {
554  // Is our absolute residual low enough?
555  if (current_residual < absolute_residual_tolerance)
556  {
557  libMesh::out << " Nonlinear solver converged, step " << step_num
558  << ", residual " << current_residual
559  << std::endl;
560  }
562  {
563  if (verbose)
564  libMesh::out << " Nonlinear solver current_residual "
565  << current_residual << " > "
566  << (absolute_residual_tolerance) << std::endl;
567  }
568 
569  // Is our relative residual low enough?
570  if ((current_residual / max_residual_norm) <
572  {
573  libMesh::out << " Nonlinear solver converged, step " << step_num
574  << ", residual reduction "
575  << current_residual / max_residual_norm
576  << " < " << relative_residual_tolerance
577  << std::endl;
578  }
580  {
581  if (verbose)
582  libMesh::out << " Nonlinear solver relative residual "
583  << (current_residual / max_residual_norm)
584  << " > " << relative_residual_tolerance
585  << std::endl;
586  }
587 
588  // For incomplete linear solves, it's not safe to test step sizes
589  if (!linear_solve_finished)
590  return;
591 
592  // Is our absolute Newton step size small enough?
593  if (step_norm < absolute_step_tolerance)
594  {
595  libMesh::out << " Nonlinear solver converged, step " << step_num
596  << ", absolute step size "
597  << step_norm
598  << " < " << absolute_step_tolerance
599  << std::endl;
600  }
601  else if (absolute_step_tolerance)
602  {
603  if (verbose)
604  libMesh::out << " Nonlinear solver absolute step size "
605  << step_norm
606  << " > " << absolute_step_tolerance
607  << std::endl;
608  }
609 
610  // Is our relative Newton step size small enough?
611  if (step_norm / max_solution_norm <
613  {
614  libMesh::out << " Nonlinear solver converged, step " << step_num
615  << ", relative step size "
616  << (step_norm / max_solution_norm)
617  << " < " << relative_step_tolerance
618  << std::endl;
619  }
620  else if (relative_step_tolerance)
621  {
622  if (verbose)
623  libMesh::out << " Nonlinear solver relative step size "
624  << (step_norm / max_solution_norm)
625  << " > " << relative_step_tolerance
626  << std::endl;
627  }
628 }
629 
630 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo