Poster of Linux kernelThe best gift for a Linux geek
NewtonSolver

NewtonSolver

Section: C Library Functions (3) Updated: Thu Apr 7 2011
Local index Up
 

NAME

NewtonSolver -  

SYNOPSIS


#include <newton_solver.h>

Inherits DiffSolver.  

Public Types


typedef DiffSolver Parent

enum SolveResult { INVALID_SOLVE_RESULT = 0, CONVERGED_NO_REASON = 1, CONVERGED_ABSOLUTE_RESIDUAL = 2, CONVERGED_RELATIVE_RESIDUAL = 4, CONVERGED_ABSOLUTE_STEP = 8, CONVERGED_RELATIVE_STEP = 16, DIVERGED_NO_REASON = 32, DIVERGED_MAX_NONLINEAR_ITERATIONS = 64, DIVERGED_BACKTRACKING_FAILURE = 128 }

typedef DifferentiableSystem sys_type
 

Public Member Functions


NewtonSolver (sys_type &system)

virtual ~NewtonSolver ()

virtual void reinit ()

virtual unsigned int solve ()

virtual unsigned int adjoint_solve ()

virtual void init ()

unsigned int total_outer_iterations ()

unsigned int total_inner_iterations ()

unsigned int solve_result ()

const sys_type & system () const

sys_type & system ()
 

Static Public Member Functions


static AutoPtr< DiffSolver > build (sys_type &s)

static std::string get_info ()

static void print_info ()

static unsigned int n_objects ()
 

Public Attributes


bool require_residual_reduction

bool brent_line_search

Real minsteplength

Real linear_tolerance_multiplier

unsigned int max_linear_iterations

unsigned int max_nonlinear_iterations

bool quiet

bool continue_after_max_iterations

bool continue_after_backtrack_failure

Real absolute_residual_tolerance

Real relative_residual_tolerance

Real absolute_step_tolerance

Real relative_step_tolerance

Real initial_linear_tolerance

Real minimum_linear_tolerance
 

Protected Types


typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 

Protected Member Functions


Real line_search (Real tol, Real last_residual, Real &current_residual, NumericVector< Number > &newton_iterate, const NumericVector< Number > &linear_solution)

void print_convergence (unsigned int step_num, Real current_residual, Real step_norm, bool linear_solve_finished)

bool test_convergence (Real current_residual, Real step_norm, bool linear_solve_finished)

void increment_constructor_count (const std::string &name)

void increment_destructor_count (const std::string &name)
 

Protected Attributes


AutoPtr< LinearSolver< Number > > linear_solver

Real max_solution_norm

Real max_residual_norm

unsigned int _outer_iterations

unsigned int _inner_iterations

sys_type & _system

unsigned int _solve_result
 

Static Protected Attributes


static Counts _counts

static Threads::atomic< unsigned int > _n_objects

static Threads::spin_mutex _mutex
 

Detailed Description

This class defines a solver which uses the default libMesh linear solver in a quasiNewton method to handle a DifferentiableSystem

This class is part of the new DifferentiableSystem framework, which is still experimental. Users of this framework should beware of bugs and future API changes.

Author:

Roy H. Stogner 2006

Definition at line 47 of file newton_solver.h.  

Member Typedef Documentation

 

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]Data structure to log the information. The log is identified by the class name.

Definition at line 105 of file reference_counter.h.  

typedef DiffSolver NewtonSolver::Parent

Definition at line 61 of file newton_solver.h.  

typedef DifferentiableSystem DiffSolver::sys_type [inherited]The type of system

Definition at line 56 of file diff_solver.h.  

Member Enumeration Documentation

 

enum DiffSolver::SolveResult [inherited]Enumeration return type for the solve() function. Multiple SolveResults may be combined (OR'd) in the single return. To test which ones are present, just AND the return value with any of the SolveResult flags defined below.

Enumerator:

INVALID_SOLVE_RESULT
A default or invalid solve result. This usually means no solve has occurred yet.
CONVERGED_NO_REASON
The solver converged but no particular reason is specified.
CONVERGED_ABSOLUTE_RESIDUAL
The DiffSolver achieved the desired absolute residual tolerance.
CONVERGED_RELATIVE_RESIDUAL
The DiffSolver achieved the desired relative residual tolerance.
CONVERGED_ABSOLUTE_STEP
The DiffSolver achieved the desired absolute step size tolerance.
CONVERGED_RELATIVE_STEP
The DiffSolver achieved the desired relative step size tolerance.
DIVERGED_NO_REASON
The DiffSolver diverged but no particular reason is specified.
DIVERGED_MAX_NONLINEAR_ITERATIONS
The DiffSolver reached the maximum allowed number of nonlinear iterations before satisfying any convergence tests.
DIVERGED_BACKTRACKING_FAILURE
The DiffSolver failed to find a descent direction by backtracking (See newton_solver.C)

Definition at line 200 of file diff_solver.h.

                   {
    INVALID_SOLVE_RESULT = 0,
    
    CONVERGED_NO_REASON = 1,

    CONVERGED_ABSOLUTE_RESIDUAL = 2,

    CONVERGED_RELATIVE_RESIDUAL = 4,

    CONVERGED_ABSOLUTE_STEP = 8,

    CONVERGED_RELATIVE_STEP = 16,

    DIVERGED_NO_REASON = 32,

    DIVERGED_MAX_NONLINEAR_ITERATIONS = 64,

    DIVERGED_BACKTRACKING_FAILURE = 128
  };
 

Constructor & Destructor Documentation

 

NewtonSolver::NewtonSolver (sys_type &system)Constructor. Requires a reference to the system to be solved.

Definition at line 228 of file newton_solver.C.

  : 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 () [virtual]Destructor.

Definition at line 240 of file newton_solver.C.

{
}
 

Member Function Documentation

 

unsigned int NewtonSolver::adjoint_solve () [virtual]This method performs an adjoint solve, using a Krylov method.

Implements DiffSolver.

Definition at line 496 of file newton_solver.C.

References DiffSolver::_system, System::add_vector(), ExplicitSystem::assemble_qoi_derivative(), DifferentiableSystem::assembly(), SparseMatrix< T >::close(), DiffSolver::CONVERGED_RELATIVE_RESIDUAL, DofMap::enforce_constraints_exactly(), System::get_dof_map(), ImplicitSystem::get_matrix(), SparseMatrix< T >::get_transpose(), ImplicitSystem::have_matrix(), linear_solver, ImplicitSystem::matrix, DiffSolver::max_linear_iterations, DiffSolver::quiet, DiffSolver::relative_residual_tolerance, ExplicitSystem::rhs, and System::update().

{
  // The adjoint_solve API (and all APIs using it) are about to see a
  // series of non-backwards-compatible changes, primarily to add
  // multiple-QoI support
  libmesh_experimental();

  START_LOG('adjoint_solve()', 'NewtonSolver');

  if (!quiet)
  std::cout << 'Assembling the Adjoint System' << std::endl;

  // Adding an adjoint_solution vector, allocate an adjoint_solution if it doesn't already exist

  NumericVector<Number> & adjoint_solution = _system.add_vector('adjoint_solution');
  
  // Do DiffSystem assembly
  _system.assembly(false, true);
  _system.matrix->close();

  // But take the adjoint
  _system.matrix->get_transpose(*_system.matrix);

  if (_system.have_matrix('Preconditioner'))
    {
      SparseMatrix<Number> &pre = _system.get_matrix('Preconditioner');
      pre.get_transpose(pre);
    }

  // And set the right hand side to the quantity of interest's
  // derivative
  _system.assemble_qoi_derivative();

  if (!quiet)
    std::cout << 'Linear solve of Adjoint System starting, tolerance ' 
              << relative_residual_tolerance << ', max iterations '
              << max_linear_iterations << std::endl;

  // Solve the linear system.  Two cases:
  const std::pair<unsigned int, Real> rval =
    (_system.have_matrix('Preconditioner')) ?
  // 1.) User-supplied preconditioner
    linear_solver->solve (*_system.matrix,
                          _system.get_matrix('Preconditioner'),
                          adjoint_solution, *_system.rhs,
                          relative_residual_tolerance,
                          max_linear_iterations) :
  // 2.) Use system matrix for the preconditioner
    linear_solver->solve (*_system.matrix,
                          adjoint_solution, *_system.rhs,
                          relative_residual_tolerance, 
                          max_linear_iterations);

  // We may need to localize a parallel solution
  _system.update ();

  // The linear solver may not have fit our constraints exactly
#ifdef LIBMESH_ENABLE_AMR
  _system.get_dof_map().enforce_constraints_exactly(_system, &adjoint_solution);
#endif

  const unsigned int linear_steps = rval.first;
  libmesh_assert(linear_steps <= max_linear_iterations);

  if (!quiet)
    std::cout << 'Linear solve of Adjoint System finished, step ' << linear_steps
              << ', residual ' << rval.second
              << std::endl;

  STOP_LOG('adjoint_solve()', 'NewtonSolver');

  // FIXME - We'll worry about getting the solve result right later...
  
  return DiffSolver::CONVERGED_RELATIVE_RESIDUAL;
}
 

AutoPtr< DiffSolver > DiffSolver::build (sys_type &s) [static, inherited]Factory. Requires a reference to the system to be solved. Returns a NewtonSolver by default

Definition at line 47 of file diff_solver.C.

Referenced by PetscDiffSolver::adjoint_solve(), and TimeSolver::init().

{
  return AutoPtr<DiffSolver>(new NewtonSolver(s));
}
 

std::string ReferenceCounter::get_info () [static, inherited]Gets a string containing the reference information.

Definition at line 45 of file reference_counter.C.

References ReferenceCounter::_counts, and Quality::name().

Referenced by ReferenceCounter::print_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)

  std::ostringstream out;
  
  out << '
      << ' ---------------------------------------------------------------------------- 
      << '| Reference count information                                                |
      << ' ---------------------------------------------------------------------------- ;
  
  for (Counts::iterator it = _counts.begin();
       it != _counts.end(); ++it)
    {
      const std::string name(it->first);
      const unsigned int creations    = it->second.first;
      const unsigned int destructions = it->second.second;

      out << '| ' << name << ' reference count information:
          << '|  Creations:    ' << creations    << '
          << '|  Destructions: ' << destructions << ';
    }
  
  out << ' ---------------------------------------------------------------------------- ;

  return out.str();

#else

  return '';
  
#endif
}
 

void ReferenceCounter::increment_constructor_count (const std::string &name) [inline, protected, inherited]Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 149 of file reference_counter.h.

References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.

Referenced by ReferenceCountedObject< Value >::ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.first++;
}
 

void ReferenceCounter::increment_destructor_count (const std::string &name) [inline, protected, inherited]Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 167 of file reference_counter.h.

References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.

Referenced by ReferenceCountedObject< Value >::~ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.second++;
}
 

void DiffSolver::init () [virtual, inherited]The initialization function. This method is used to initialize internal data structures before a simulation begins.

Reimplemented in PetscDiffSolver.

Definition at line 63 of file diff_solver.C.

References DiffSolver::max_residual_norm, and DiffSolver::max_solution_norm.

Referenced by PetscDiffSolver::init().

{
  // Reset the max_step_size and max_residual_norm for a new problem
  max_solution_norm = 0.;
  max_residual_norm = 0.;
}
 

Real NewtonSolver::line_search (Realtol, Reallast_residual, Real &current_residual, NumericVector< Number > &newton_iterate, const NumericVector< Number > &linear_solution) [protected]This does a line search in the direction opposite linear_solution to try and minimize the residual of newton_iterate. newton_iterate is moved to the end of the quasiNewton step, and the return value is the substep size.

Definition at line 40 of file newton_solver.C.

References DiffSolver::_outer_iterations, DiffSolver::_solve_result, DiffSolver::_system, NumericVector< T >::add(), DifferentiableSystem::assembly(), brent_line_search, NumericVector< T >::close(), DiffSolver::continue_after_backtrack_failure, DiffSolver::DIVERGED_BACKTRACKING_FAILURE, NumericVector< T >::l2_norm(), std::max(), std::min(), minsteplength, DiffSolver::quiet, require_residual_reduction, ExplicitSystem::rhs, and SIGN().

Referenced by solve().

{
  // 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_convergence_failure();
            }
          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 = current_residual = 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;
}
 

static unsigned int ReferenceCounter::n_objects () [inline, static, inherited]Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 76 of file reference_counter.h.

References ReferenceCounter::_n_objects.

Referenced by System::read_serialized_blocked_dof_objects(), and System::write_serialized_blocked_dof_objects().

  { return _n_objects; }
 

void NewtonSolver::print_convergence (unsigned intstep_num, Realcurrent_residual, Realstep_norm, boollinear_solve_finished) [protected]This prints output for the convergence criteria based on by the given residual and step size.

Definition at line 621 of file newton_solver.C.

References DiffSolver::absolute_residual_tolerance, DiffSolver::absolute_step_tolerance, DiffSolver::max_residual_norm, DiffSolver::max_solution_norm, DiffSolver::quiet, DiffSolver::relative_residual_tolerance, and DiffSolver::relative_step_tolerance.

Referenced by solve().

{
  // Is our absolute residual low enough?
  if (current_residual < absolute_residual_tolerance)
    {
      std::cout << '  Nonlinear solver converged, step ' << step_num
                << ', residual ' << current_residual
                << std::endl;
    }
  else if (absolute_residual_tolerance)
    {
      std::cout << '  Nonlinear solver current_residual '
                << current_residual << ' > '
                << (absolute_residual_tolerance) << std::endl;
    }

  // Is our relative residual low enough?
  if ((current_residual / max_residual_norm) <
      relative_residual_tolerance)
    {
      std::cout << '  Nonlinear solver converged, step ' << step_num
                << ', residual reduction '
                << current_residual / max_residual_norm
                << ' < ' << relative_residual_tolerance
                << std::endl;
    }
  else if (relative_residual_tolerance)
    {
      if (!quiet)
        std::cout << '  Nonlinear solver relative residual '
                  << (current_residual / max_residual_norm)
                  << ' > ' << relative_residual_tolerance
                  << std::endl;
    }

  // For incomplete linear solves, it's not safe to test step sizes
  if (!linear_solve_finished)
    return;

  // Is our absolute Newton step size small enough?
  if (step_norm < absolute_step_tolerance)
    {
      std::cout << '  Nonlinear solver converged, step ' << step_num
                << ', absolute step size '
                << step_norm
                << ' < ' << absolute_step_tolerance
                << std::endl;
    }
  else if (absolute_step_tolerance)
    {
      std::cout << '  Nonlinear solver absolute step size '
                << step_norm
                << ' > ' << absolute_step_tolerance
                << std::endl;
    }

  // Is our relative Newton step size small enough?
  if (step_norm / max_solution_norm <
      relative_step_tolerance)
    {
      std::cout << '  Nonlinear solver converged, step ' << step_num
                << ', relative step size '
                << (step_norm / max_solution_norm)
                << ' < ' << relative_step_tolerance
                << std::endl;
    }
  else if (relative_step_tolerance)
    {
      std::cout << '  Nonlinear solver relative step size '
                << (step_norm / max_solution_norm)
                << ' > ' << relative_step_tolerance
                << std::endl;
    }
}
 

void ReferenceCounter::print_info () [static, inherited]Prints the reference information to std::cout.

Definition at line 83 of file reference_counter.C.

References ReferenceCounter::get_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
  
  std::cout << ReferenceCounter::get_info();
  
#endif
}
 

void NewtonSolver::reinit () [virtual]The reinitialization function. This method is used after changes in the mesh.

Reimplemented from DiffSolver.

Definition at line 246 of file newton_solver.C.

References linear_solver, and DiffSolver::reinit().

{
  Parent::reinit();

  linear_solver->clear();
}
 

unsigned int NewtonSolver::solve () [virtual]This method performs a solve, using an inexact Newton-Krylov method with line search.

Implements DiffSolver.

Definition at line 255 of file newton_solver.C.

References DiffSolver::_inner_iterations, DiffSolver::_outer_iterations, DiffSolver::_solve_result, DiffSolver::_system, NumericVector< T >::add(), DifferentiableSystem::assembly(), NumericVector< T >::clone(), NumericVector< T >::close(), DiffSolver::continue_after_max_iterations, DiffSolver::DIVERGED_BACKTRACKING_FAILURE, DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS, DofMap::enforce_constraints_exactly(), System::get_dof_map(), ImplicitSystem::get_matrix(), ImplicitSystem::have_matrix(), DiffSolver::initial_linear_tolerance, DiffSolver::INVALID_SOLVE_RESULT, NumericVector< T >::l2_norm(), line_search(), linear_solver, linear_tolerance_multiplier, ImplicitSystem::matrix, std::max(), DiffSolver::max_linear_iterations, DiffSolver::max_nonlinear_iterations, DiffSolver::max_residual_norm, DiffSolver::max_solution_norm, std::min(), DiffSolver::minimum_linear_tolerance, print_convergence(), DiffSolver::quiet, ExplicitSystem::rhs, System::solution, test_convergence(), System::update(), and NumericVector< T >::zero().

{
  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 LIBMESH_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 the 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_convergence_failure();
          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: '
                  << current_residual << std::endl;

      // Make sure our linear tolerance is low enough
      current_linear_tolerance = std::min (current_linear_tolerance,
        current_residual * linear_tolerance_multiplier);

      // But don't let it be too small
      if (current_linear_tolerance < minimum_linear_tolerance)
        {
          current_linear_tolerance = minimum_linear_tolerance;
        }

      // At this point newton_iterate is the current guess, and
      // linear_solution is now about to become the NEGATIVE of the next
      // Newton step.

      // Our best initial guess for the linear_solution is zero!
      linear_solution.zero();

      if (!quiet)
        std::cout << 'Linear solve starting, tolerance ' 
                  << current_linear_tolerance << std::endl;

      // Solve the linear system.  Two cases:
      const std::pair<unsigned int, Real> rval =
        (_system.have_matrix('Preconditioner')) ?
      // 1.) User-supplied preconditioner
        linear_solver->solve (matrix, _system.get_matrix('Preconditioner'),
                              linear_solution, rhs, current_linear_tolerance,
                              max_linear_iterations) :
      // 2.) Use system matrix for the preconditioner
        linear_solver->solve (matrix, linear_solution, rhs,
                              current_linear_tolerance, 
                              max_linear_iterations);
      // We may need to localize a parallel solution
      _system.update ();
      // The linear solver may not have fit our constraints exactly
#ifdef LIBMESH_ENABLE_AMR
      _system.get_dof_map().enforce_constraints_exactly(_system, &linear_solution);
#endif

      const unsigned int linear_steps = rval.first;
      libmesh_assert(linear_steps <= max_linear_iterations);
      _inner_iterations += linear_steps;

      const bool linear_solve_finished = 
        !(linear_steps == max_linear_iterations);

      if (!quiet)
        std::cout << 'Linear solve finished, step ' << linear_steps
                  << ', residual ' << rval.second
                  << std::endl;

      // Compute the l2 norm of the nonlinear update
      Real norm_delta = linear_solution.l2_norm();

      if (!quiet)
        std::cout << 'Trying full Newton step' << std::endl;
      // Take a full Newton step
      newton_iterate.add (-1., linear_solution);
      newton_iterate.close();

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

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

// A potential method for avoiding oversolving?
/*
      Real predicted_absolute_error =
        current_residual * norm_delta / last_residual;

      Real predicted_relative_error =
        predicted_absolute_error / max_solution_norm;

      std::cout << 'Predicted absolute error = ' <<
        predicted_absolute_error << std::endl;

      std::cout << 'Predicted relative error = ' <<
        predicted_relative_error << std::endl;
*/

      // don't fiddle around if we've already converged
      if (test_convergence(current_residual, norm_delta,
                           linear_solve_finished))
        {
          if (!quiet)
            print_convergence(_outer_iterations, current_residual,
                              norm_delta, linear_solve_finished);
          _outer_iterations++;
          break; // out of _outer_iterations for loop
        }

      // otherwise, backtrack if necessary
      Real steplength =
        this->line_search(std::sqrt(TOLERANCE),
                          last_residual, current_residual,
                          newton_iterate, linear_solution);
      norm_delta *= steplength;

      // Check to see if backtracking failed,
      // and break out of the nonlinear loop if so...
      if (_solve_result == DiffSolver::DIVERGED_BACKTRACKING_FAILURE)
        {
          _outer_iterations++;
          break; // out of _outer_iterations for loop
        }

      // Compute the l2 norm of the whole solution
      norm_total = newton_iterate.l2_norm();

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

      // Print out information for the 
      // nonlinear iterations.
      if (!quiet)
        std::cout << '  Nonlinear step: |du|/|u| = '
                  << norm_delta / norm_total
                  << ', |du| = ' << norm_delta
                  << std::endl;

      // Terminate the solution iteration if the difference between
      // this iteration and the last is sufficiently small.
      if (!quiet)
        print_convergence(_outer_iterations, current_residual,
                          norm_delta / steplength,
                          linear_solve_finished);
      if (test_convergence(current_residual, norm_delta / steplength,
                           linear_solve_finished))
        {
          _outer_iterations++;
          break; // out of _outer_iterations for loop
        }

      if (_outer_iterations >= max_nonlinear_iterations - 1)
        {
          std::cout << '  Nonlinear solver FAILED TO CONVERGE by step '
                    << _outer_iterations << ' with norm '
                    << norm_total << std::endl;
          if (continue_after_max_iterations)
            {
              _solve_result = DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
              std::cout << '  Continuing anyway...' << std::endl;
            }
          else
            {
              libmesh_convergence_failure();
            }
          continue;
        }
    } // end nonlinear loop

  // The linear solver may not have fit our constraints exactly
#ifdef LIBMESH_ENABLE_AMR
  _system.get_dof_map().enforce_constraints_exactly(_system);
#endif

  // We may need to localize a parallel solution
  _system.update ();

  STOP_LOG('solve()', 'NewtonSolver');

  // Make sure we are returning something sensible as the
  // _solve_result.
  libmesh_assert (_solve_result != DiffSolver::INVALID_SOLVE_RESULT);
  
  return _solve_result;
}
 

unsigned int DiffSolver::solve_result () [inline, inherited]Returns:

the value of the SolveResult from the last solve.

Definition at line 116 of file diff_solver.h.

References DiffSolver::_solve_result.

{ return _solve_result; }
 

sys_type& DiffSolver::system () [inline, inherited]Returns:

a writeable reference to the system we are solving.

Definition at line 126 of file diff_solver.h.

References DiffSolver::_system.

{ return _system; }
 

const sys_type& DiffSolver::system () const [inline, inherited]Returns:

a constant reference to the system we are solving.

Definition at line 121 of file diff_solver.h.

References DiffSolver::_system.

Referenced by __libmesh_petsc_diff_solver_jacobian(), and __libmesh_petsc_diff_solver_residual().

{ return _system; }
 

bool NewtonSolver::test_convergence (Realcurrent_residual, Realstep_norm, boollinear_solve_finished) [protected]This returns true if a convergence criterion has been passed by the given residual and step size; false otherwise.

Definition at line 574 of file newton_solver.C.

References DiffSolver::_solve_result, DiffSolver::absolute_residual_tolerance, DiffSolver::absolute_step_tolerance, DiffSolver::CONVERGED_ABSOLUTE_RESIDUAL, DiffSolver::CONVERGED_ABSOLUTE_STEP, DiffSolver::CONVERGED_RELATIVE_RESIDUAL, DiffSolver::CONVERGED_RELATIVE_STEP, DiffSolver::max_residual_norm, DiffSolver::max_solution_norm, DiffSolver::relative_residual_tolerance, and DiffSolver::relative_step_tolerance.

Referenced by solve().

{
  // We haven't converged unless we pass a convergence test
  bool has_converged = false;

  // Is our absolute residual low enough?
  if (current_residual < absolute_residual_tolerance)
    {
      _solve_result |= CONVERGED_ABSOLUTE_RESIDUAL;
      has_converged = true;
    }
  
  // Is our relative residual low enough?
  if ((current_residual / max_residual_norm) <
      relative_residual_tolerance)
    {
      _solve_result |= CONVERGED_RELATIVE_RESIDUAL;
      has_converged = true;
    }
  
  // For incomplete linear solves, it's not safe to test step sizes
  if (!linear_solve_finished)
    {
      return has_converged;
    }
  
  // Is our absolute Newton step size small enough?
  if (step_norm < absolute_step_tolerance)
    {
      _solve_result |= CONVERGED_ABSOLUTE_STEP;
      has_converged = true;
    }
  
  // Is our relative Newton step size small enough?
  if (step_norm / max_solution_norm <
      relative_step_tolerance)
    {
      _solve_result |= CONVERGED_RELATIVE_STEP;
      has_converged = true;
    }
  
  return has_converged;
}
 

unsigned int DiffSolver::total_inner_iterations () [inline, inherited]Returns:

the number of 'inner' (e.g. Krylov) iterations required by the last solve.

Definition at line 111 of file diff_solver.h.

References DiffSolver::_inner_iterations.

{ return _inner_iterations; }
 

unsigned int DiffSolver::total_outer_iterations () [inline, inherited]Returns:

the number of 'outer' (e.g. quasi-Newton) iterations required by the last solve.

Definition at line 105 of file diff_solver.h.

References DiffSolver::_outer_iterations.

{ return _outer_iterations; }
 

Member Data Documentation

 

ReferenceCounter::Counts ReferenceCounter::_counts [static, protected, inherited]Actually holds the data.

Definition at line 110 of file reference_counter.h.

Referenced by ReferenceCounter::get_info(), ReferenceCounter::increment_constructor_count(), and ReferenceCounter::increment_destructor_count().  

unsigned int DiffSolver::_inner_iterations [protected, inherited]The number of inner iterations used by the last solve

Definition at line 280 of file diff_solver.h.

Referenced by solve(), and DiffSolver::total_inner_iterations().  

Threads::spin_mutex ReferenceCounter::_mutex [static, protected, inherited]Mutual exclusion object to enable thread-safe reference counting.

Definition at line 123 of file reference_counter.h.  

Threads::atomic< unsigned int > ReferenceCounter::_n_objects [static, protected, inherited]The number of objects. Print the reference count information when the number returns to 0.

Definition at line 118 of file reference_counter.h.

Referenced by ReferenceCounter::n_objects(), ReferenceCounter::ReferenceCounter(), and ReferenceCounter::~ReferenceCounter().  

unsigned int DiffSolver::_outer_iterations [protected, inherited]The number of outer iterations used by the last solve

Definition at line 275 of file diff_solver.h.

Referenced by line_search(), PetscDiffSolver::solve(), solve(), and DiffSolver::total_outer_iterations().  

unsigned int DiffSolver::_solve_result [protected, inherited]Initialized to zero. solve_result is typically set internally in the solve() function before it returns. When non-zero, solve_result tells the result of the latest solve. See enum definition for description.

Definition at line 293 of file diff_solver.h.

Referenced by line_search(), solve(), DiffSolver::solve_result(), and test_convergence().  

sys_type& DiffSolver::_system [protected, inherited]A reference to the system we are solving.

Definition at line 285 of file diff_solver.h.

Referenced by PetscDiffSolver::adjoint_solve(), adjoint_solve(), line_search(), PetscDiffSolver::solve(), solve(), and DiffSolver::system().  

Real DiffSolver::absolute_residual_tolerance [inherited]The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolerance or less than relative_residual_tolerance times the initial residual.

Users should increase any of these tolerances that they want to use for a stopping condition.

Definition at line 169 of file diff_solver.h.

Referenced by print_convergence(), and test_convergence().  

Real DiffSolver::absolute_step_tolerance [inherited]The DiffSolver should exit after the full nonlinear step norm is reduced to either less than absolute_step_tolerance or less than relative_step_tolerance times the largest nonlinear solution which has been seen so far.

Users should increase any of these tolerances that they want to use for a stopping condition.

Definition at line 181 of file diff_solver.h.

Referenced by print_convergence(), and test_convergence().  

bool NewtonSolver::brent_line_searchIf require_residual_reduction is true, the solver may reduce step lengths when required. If so, brent_line_search is an option. If brent_line_search is set to false, the solver reduces the length of its steps by 1/2 iteratively until it finds residual reduction. If true, step lengths are first reduced by 1/2 or more to find some residual reduction, then Brent's method is used to find as much residual reduction as possible.

brent_line_search is currently set to true by default.

Definition at line 100 of file newton_solver.h.

Referenced by line_search().  

bool DiffSolver::continue_after_backtrack_failure [inherited]Defaults to false, telling the DiffSolver to throw a libmesh_error() when the backtracking scheme fails to find a descent direction.

Definition at line 158 of file diff_solver.h.

Referenced by line_search().  

bool DiffSolver::continue_after_max_iterations [inherited]Defaults to true, telling the DiffSolver to continue rather than exit when a solve has reached its maximum number of nonlinear iterations.

Definition at line 152 of file diff_solver.h.

Referenced by solve().  

Real DiffSolver::initial_linear_tolerance [inherited]Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the tolerance for later solves.

Definition at line 188 of file diff_solver.h.

Referenced by solve().  

AutoPtr<LinearSolver<Number> > NewtonSolver::linear_solver [protected]The LinearSolver defines the interface used to solve the linear_implicit system. This class handles all the details of interfacing with various linear algebra packages like PETSc or LASPACK.

Definition at line 124 of file newton_solver.h.

Referenced by adjoint_solve(), reinit(), and solve().  

Real NewtonSolver::linear_tolerance_multiplierThe tolerance for linear solves is kept below this multiplier (which defaults to 1e-3) times the norm of the current nonlinear residual

Definition at line 114 of file newton_solver.h.

Referenced by solve().  

unsigned int DiffSolver::max_linear_iterations [inherited]Each linear solver step should exit after max_linear_iterations is exceeded.

Definition at line 132 of file diff_solver.h.

Referenced by PetscDiffSolver::adjoint_solve(), adjoint_solve(), ContinuationSystem::continuation_solve(), solve(), and ContinuationSystem::solve_tangent().  

unsigned int DiffSolver::max_nonlinear_iterations [inherited]The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_iterations is false, or should end the nonlinear solve if max_nonlinear_iterations is exceeded and continue_after_max_iterations is true.

Definition at line 140 of file diff_solver.h.

Referenced by ContinuationSystem::continuation_solve(), solve(), and ContinuationSystem::update_solution().  

Real DiffSolver::max_residual_norm [protected, inherited]The largest nonlinear residual which the DiffSolver has yet seen will be stored here, to be used for stopping criteria based on relative_residual_tolerance

Definition at line 270 of file diff_solver.h.

Referenced by DiffSolver::init(), print_convergence(), DiffSolver::reinit(), solve(), and test_convergence().  

Real DiffSolver::max_solution_norm [protected, inherited]The largest solution norm which the DiffSolver has yet seen will be stored here, to be used for stopping criteria based on relative_step_tolerance

Definition at line 263 of file diff_solver.h.

Referenced by DiffSolver::init(), print_convergence(), DiffSolver::reinit(), solve(), and test_convergence().  

Real DiffSolver::minimum_linear_tolerance [inherited]The tolerance for linear solves is kept above this minimum

Definition at line 193 of file diff_solver.h.

Referenced by solve().  

Real NewtonSolver::minsteplengthIf the quasi-Newton step length must be reduced to below this factor to give a residual reduction, then the Newton solver dies with a libmesh_error() It is currently set to 1e-5 by default.

Definition at line 108 of file newton_solver.h.

Referenced by line_search().  

bool DiffSolver::quiet [inherited]The DiffSolver should not print anything to std::cout unless quiet is set to false

Definition at line 146 of file diff_solver.h.

Referenced by adjoint_solve(), line_search(), print_convergence(), and solve().  

Real DiffSolver::relative_residual_tolerance [inherited]

Definition at line 170 of file diff_solver.h.

Referenced by PetscDiffSolver::adjoint_solve(), adjoint_solve(), print_convergence(), and test_convergence().  

Real DiffSolver::relative_step_tolerance [inherited]

Definition at line 182 of file diff_solver.h.

Referenced by print_convergence(), and test_convergence().  

bool NewtonSolver::require_residual_reductionIf this is set to true, the solver is forced to test the residual after each Newton step, and to reduce the length of its steps whenever necessary to avoid a residual increase. It is currently set to true by default; set it to false to avoid unnecessary residual assembly on well-behaved systems.

Definition at line 87 of file newton_solver.h.

Referenced by line_search().

 

Author

Generated automatically by Doxygen for libMesh from the source code.


 

Index

NAME
SYNOPSIS
Public Types
Public Member Functions
Static Public Member Functions
Public Attributes
Protected Types
Protected Member Functions
Protected Attributes
Static Protected Attributes
Detailed Description
Member Typedef Documentation
typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]Data structure to log the information. The log is identified by the class name.
typedef DiffSolver NewtonSolver::Parent
typedef DifferentiableSystem DiffSolver::sys_type [inherited]The type of system
Member Enumeration Documentation
enum DiffSolver::SolveResult [inherited]Enumeration return type for the solve() function. Multiple SolveResults may be combined (OR'd) in the single return. To test which ones are present, just AND the return value with any of the SolveResult flags defined below.
Constructor & Destructor Documentation
NewtonSolver::NewtonSolver (sys_type &system)Constructor. Requires a reference to the system to be solved.
NewtonSolver::~NewtonSolver () [virtual]Destructor.
Member Function Documentation
unsigned int NewtonSolver::adjoint_solve () [virtual]This method performs an adjoint solve, using a Krylov method.
AutoPtr< DiffSolver > DiffSolver::build (sys_type &s) [static, inherited]Factory. Requires a reference to the system to be solved. Returns a NewtonSolver by default
std::string ReferenceCounter::get_info () [static, inherited]Gets a string containing the reference information.
void ReferenceCounter::increment_constructor_count (const std::string &name) [inline, protected, inherited]Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.
void ReferenceCounter::increment_destructor_count (const std::string &name) [inline, protected, inherited]Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.
void DiffSolver::init () [virtual, inherited]The initialization function. This method is used to initialize internal data structures before a simulation begins.
Real NewtonSolver::line_search (Realtol, Reallast_residual, Real &current_residual, NumericVector< Number > &newton_iterate, const NumericVector< Number > &linear_solution) [protected]This does a line search in the direction opposite linear_solution to try and minimize the residual of newton_iterate. newton_iterate is moved to the end of the quasiNewton step, and the return value is the substep size.
static unsigned int ReferenceCounter::n_objects () [inline, static, inherited]Prints the number of outstanding (created, but not yet destroyed) objects.
void NewtonSolver::print_convergence (unsigned intstep_num, Realcurrent_residual, Realstep_norm, boollinear_solve_finished) [protected]This prints output for the convergence criteria based on by the given residual and step size.
void ReferenceCounter::print_info () [static, inherited]Prints the reference information to std::cout.
void NewtonSolver::reinit () [virtual]The reinitialization function. This method is used after changes in the mesh.
unsigned int NewtonSolver::solve () [virtual]This method performs a solve, using an inexact Newton-Krylov method with line search.
unsigned int DiffSolver::solve_result () [inline, inherited]Returns:
sys_type& DiffSolver::system () [inline, inherited]Returns:
const sys_type& DiffSolver::system () const [inline, inherited]Returns:
bool NewtonSolver::test_convergence (Realcurrent_residual, Realstep_norm, boollinear_solve_finished) [protected]This returns true if a convergence criterion has been passed by the given residual and step size; false otherwise.
unsigned int DiffSolver::total_inner_iterations () [inline, inherited]Returns:
unsigned int DiffSolver::total_outer_iterations () [inline, inherited]Returns:
Member Data Documentation
ReferenceCounter::Counts ReferenceCounter::_counts [static, protected, inherited]Actually holds the data.
unsigned int DiffSolver::_inner_iterations [protected, inherited]The number of inner iterations used by the last solve
Threads::spin_mutex ReferenceCounter::_mutex [static, protected, inherited]Mutual exclusion object to enable thread-safe reference counting.
Threads::atomic< unsigned int > ReferenceCounter::_n_objects [static, protected, inherited]The number of objects. Print the reference count information when the number returns to 0.
unsigned int DiffSolver::_outer_iterations [protected, inherited]The number of outer iterations used by the last solve
unsigned int DiffSolver::_solve_result [protected, inherited]Initialized to zero. solve_result is typically set internally in the solve() function before it returns. When non-zero, solve_result tells the result of the latest solve. See enum definition for description.
sys_type& DiffSolver::_system [protected, inherited]A reference to the system we are solving.
Real DiffSolver::absolute_residual_tolerance [inherited]The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolerance or less than relative_residual_tolerance times the initial residual.
Real DiffSolver::absolute_step_tolerance [inherited]The DiffSolver should exit after the full nonlinear step norm is reduced to either less than absolute_step_tolerance or less than relative_step_tolerance times the largest nonlinear solution which has been seen so far.
bool NewtonSolver::brent_line_searchIf require_residual_reduction is true, the solver may reduce step lengths when required. If so, brent_line_search is an option. If brent_line_search is set to false, the solver reduces the length of its steps by 1/2 iteratively until it finds residual reduction. If true, step lengths are first reduced by 1/2 or more to find some residual reduction, then Brent's method is used to find as much residual reduction as possible.
bool DiffSolver::continue_after_backtrack_failure [inherited]Defaults to false, telling the DiffSolver to throw a libmesh_error() when the backtracking scheme fails to find a descent direction.
bool DiffSolver::continue_after_max_iterations [inherited]Defaults to true, telling the DiffSolver to continue rather than exit when a solve has reached its maximum number of nonlinear iterations.
Real DiffSolver::initial_linear_tolerance [inherited]Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the tolerance for later solves.
AutoPtr<LinearSolver<Number> > NewtonSolver::linear_solver [protected]The LinearSolver defines the interface used to solve the linear_implicit system. This class handles all the details of interfacing with various linear algebra packages like PETSc or LASPACK.
Real NewtonSolver::linear_tolerance_multiplierThe tolerance for linear solves is kept below this multiplier (which defaults to 1e-3) times the norm of the current nonlinear residual
unsigned int DiffSolver::max_linear_iterations [inherited]Each linear solver step should exit after max_linear_iterations is exceeded.
unsigned int DiffSolver::max_nonlinear_iterations [inherited]The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_iterations is false, or should end the nonlinear solve if max_nonlinear_iterations is exceeded and continue_after_max_iterations is true.
Real DiffSolver::max_residual_norm [protected, inherited]The largest nonlinear residual which the DiffSolver has yet seen will be stored here, to be used for stopping criteria based on relative_residual_tolerance
Real DiffSolver::max_solution_norm [protected, inherited]The largest solution norm which the DiffSolver has yet seen will be stored here, to be used for stopping criteria based on relative_step_tolerance
Real DiffSolver::minimum_linear_tolerance [inherited]The tolerance for linear solves is kept above this minimum
Real NewtonSolver::minsteplengthIf the quasi-Newton step length must be reduced to below this factor to give a residual reduction, then the Newton solver dies with a libmesh_error() It is currently set to 1e-5 by default.
bool DiffSolver::quiet [inherited]The DiffSolver should not print anything to std::cout unless quiet is set to false
Real DiffSolver::relative_residual_tolerance [inherited]
Real DiffSolver::relative_step_tolerance [inherited]
bool NewtonSolver::require_residual_reductionIf this is set to true, the solver is forced to test the residual after each Newton step, and to reduce the length of its steps whenever necessary to avoid a residual increase. It is currently set to true by default; set it to false to avoid unnecessary residual assembly on well-behaved systems.
Author

This document was created by man2html, using the manual pages.
Time: 21:51:24 GMT, April 16, 2011