newton_solver.c

一个用来实现偏微分方程中网格的计算库

#include "cmath" // For isnan(), when it's defined

#include "diff_system.h"
#include "equation_systems.h"
#include "libmesh_logging.h"
#include "linear_solver.h"
#include "newton_solver.h"
#include "numeric_vector.h"
#include "sparse_matrix.h"
#include "dof_map.h"

// SIGN from Numerical Recipes
template <typename T>
inline
T SIGN(T a, T b)
{
  return b >= 0 ? std::abs(a) : -std::abs(a);
}

Real NewtonSolver::line_search(Real tol,
                               Real last_residual,
                               Real current_residual,
                               NumericVector<Number> &newton_iterate,
                               const NumericVector<Number> &linear_solution)
{
  // Take a full step if we got a residual reduction or if we
  // aren't substepping
  if (current_residual < last_residual ||
      !require_residual_reduction)
    return 1.;

  // The residual vector
  NumericVector<Number> &rhs = *(_system.rhs);

  Real ax = 0.;  // First abscissa, don't take negative steps
  Real cx = 1.;  // Second abscissa, don't extrapolate steps

  // Find bx, a step length that gives lower residual than ax or cx
  Real bx = 1.;
  while (current_residual > last_residual)
    {
      // Reduce step size to 1/2, 1/4, etc.
      Real substepdivision;
      if (brent_line_search)
        {
          substepdivision = std::min(0.5, last_residual/current_residual);
          substepdivision = std::max(substepdivision, tol*2.);
        }
      else
        substepdivision = 0.5;

      newton_iterate.add (bx * (1.-substepdivision),
                          linear_solution);
      newton_iterate.close();
      bx *= substepdivision;
      if (!quiet)
        std::cout << "  Shrinking Newton step to "
                  << bx << std::endl;

      // Check residual with fractional Newton step
      _system.assembly (true, false);

      rhs.close();
      current_residual = rhs.l2_norm();
      if (!quiet)
        std::cout << "  Current Residual: "
                  << current_residual << std::endl;

      if (bx/2. < minsteplength && 
          current_residual > last_residual)
        {
          std::cout << "Inexact Newton step FAILED at step "
                    << _outer_iterations << std::endl;

          if (!continue_after_backtrack_failure)
            {
              libmesh_error();
            }
          else
            {
              std::cout << "Continuing anyway ..." << std::endl;
              _solve_result = DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
              return bx;
            }
        }
    } // end while (current_residual > last_residual)

  // Now return that reduced-residual step, or  use Brent's method to
  // find a more optimal step.

  if (!brent_line_search)
    return bx;

  // Brent's method adapted from Numerical Recipes in C, ch. 10.2
  Real e = 0.;

  Real x = bx, w = bx, v = bx;

  // Residuals at bx
  Real fx = current_residual,
       fw = current_residual,
       fv = current_residual;

  // Max iterations for Brent's method loop
  const unsigned int max_i = 20;

  // for golden ratio steps
  const Real golden_ratio = 1.-(std::sqrt(5.)-1.)/2.;

  for (unsigned int i=1; i <= max_i; i++)
    {
      Real xm = (ax+cx)*0.5;
      Real tol1 = tol * std::abs(x) + tol*tol;
      Real tol2 = 2.0 * tol1;

      // Test if we're done
      if (std::abs(x-xm) <= (tol2 - 0.5 * (cx - ax)))
        return x;

      Real d;

      // Construct a parabolic fit
      if (std::abs(e) > tol1)
        {
          Real r = (x-w)*(fx-fv);
          Real q = (x-v)*(fx-fw);
          Real p = (x-v)*q-(x-w)*r;
          q = 2. * (q-r);
          if (q > 0.)
            p = -p;
          else
            q = std::abs(q);
          if (std::abs(p) >= std::abs(0.5*q*e) ||
              p <= q * (ax-x) ||
              p >= q * (cx-x))
            {
              // Take a golden section step
              e = x >= xm ? ax-x : cx-x;
              d = golden_ratio * e;
            }
          else
            {
              // Take a parabolic fit step
              d = p/q;
              if (x+d-ax < tol2 || cx-(x+d) < tol2)
                d = SIGN(tol1, xm - x);
            }
        }
      else
        {
          // Take a golden section step
          e = x >= xm ? ax-x : cx-x;
          d = golden_ratio * e;
        }

      Real u = std::abs(d) >= tol1 ? x+d : x + SIGN(tol1,d);

      // Assemble the residual at the new steplength u
      newton_iterate.add (bx - u, linear_solution);
      newton_iterate.close();
      bx = u;
      if (!quiet)
        std::cout << "  Shrinking Newton step to "
                  << bx << std::endl;

      _system.assembly (true, false);

      rhs.close();
      Real fu = rhs.l2_norm();
      if (!quiet)
        std::cout << "  Current Residual: "
                  << fu << std::endl;

      if (fu <= fx)
        {
          if (u >= x)
            ax = x;
          else
            cx = x;
          v = w;   w = x;   x = u;
          fv = fw; fw = fx; fx = fu;
        }
      else
        {
          if (u < x)
            ax = u;
          else
            cx = u;
          if (fu <= fw || w == x)
            {
              v = w;   w = u;
              fv = fw; fw = fu;
            }
          else if (fu <= fv || v == x || v == w)
            {
              v = u;
              fv = fu;
            }
        }
    }

  std::cout << "Warning!  Too many iterations used in Brent line search!"
            << std::endl;
  return bx;
}


NewtonSolver::NewtonSolver (sys_type& s)
  : Parent(s),
    require_residual_reduction(true),
    brent_line_search(true),
    minsteplength(1e-5),
    linear_tolerance_multiplier(1e-3),
    linear_solver(LinearSolver<Number>::build())
{
}



NewtonSolver::~NewtonSolver ()
{
}



void NewtonSolver::reinit()
{
  Parent::reinit();

  linear_solver->clear();
}



unsigned int NewtonSolver::solve()
{
  START_LOG("solve()", "NewtonSolver");

  // Reset any prior solve result
  _solve_result = INVALID_SOLVE_RESULT;
  
  NumericVector<Number> &newton_iterate = *(_system.solution);

  AutoPtr<NumericVector<Number> > linear_solution_ptr = newton_iterate.clone();
  NumericVector<Number> &linear_solution = *linear_solution_ptr;
  NumericVector<Number> &rhs = *(_system.rhs);

  newton_iterate.close();
  linear_solution.close();
  rhs.close();

#ifdef ENABLE_AMR
  _system.get_dof_map().enforce_constraints_exactly(_system);
#endif

  SparseMatrix<Number> &matrix = *(_system.matrix);

  // Prepare to take incomplete steps
  Real last_residual=0.;

  // Set starting linear tolerance
  Real current_linear_tolerance = initial_linear_tolerance;

  // Start counting our linear solver steps
  _inner_iterations = 0;

  // Now we begin the nonlinear loop
  for (_outer_iterations=0; _outer_iterations<max_nonlinear_iterations;
       ++_outer_iterations)
    {
      if (!quiet)
        std::cout << "Assembling System" << std::endl;

      _system.assembly(true, true);
      rhs.close();
      Real current_residual = rhs.l2_norm();
      last_residual = current_residual;

// isnan() isn't standard C++ yet
#ifdef isnan
      if (isnan(current_residual))
        {
          std::cout << "  Nonlinear solver DIVERGED at step " << last_residual
                    << " with norm Not-a-Number"
                    << std::endl;
          libmesh_error();
          continue;
        }
#endif // isnan

      max_residual_norm = std::max (current_residual,
                                    max_residual_norm);
 
      // Compute the l2 norm of the whole solution
      Real norm_total = newton_iterate.l2_norm();

      max_solution_norm = std::max(max_solution_norm, norm_total);

      if (!quiet)
        std::cout << "Nonlinear Residual: "