libMesh::NewtonSolver Class Reference

#include <newton_solver.h>

Inheritance diagram for libMesh::NewtonSolver:

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, DIVERGED_LINEAR_SOLVER_FAILURE = 256
}
 
typedef ImplicitSystem sys_type
 

Public Member Functions

 NewtonSolver (sys_type &system)
 
virtual ~NewtonSolver ()
 
virtual void init () libmesh_override
 
virtual void reinit () libmesh_override
 
virtual unsigned int solve () libmesh_override
 
LinearSolver< Number > & get_linear_solver ()
 
const LinearSolver< Number > & get_linear_solver () const
 
unsigned int total_outer_iterations ()
 
unsigned int total_inner_iterations ()
 
unsigned int solve_result ()
 
const sys_typesystem () const
 
sys_typesystem ()
 
const Parallel::Communicatorcomm () const
 
processor_id_type n_processors () const
 
processor_id_type processor_id () const
 

Static Public Member Functions

static std::unique_ptr< DiffSolverbuild (sys_type &s)
 
static std::string get_info ()
 
static void print_info (std::ostream &out=libMesh::out)
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

Public Attributes

bool require_residual_reduction
 
bool require_finite_residual
 
bool brent_line_search
 
bool track_linear_convergence
 
Real minsteplength
 
Real linear_tolerance_multiplier
 
unsigned int max_linear_iterations
 
unsigned int max_nonlinear_iterations
 
bool quiet
 
bool verbose
 
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
 
std::unique_ptr< LinearSolutionMonitorlinear_solution_monitor
 

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

std::unique_ptr< 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
 
const Parallel::Communicator_communicator
 

Static Protected Attributes

static Counts _counts
 
static Threads::atomic< unsigned int > _n_objects
 
static Threads::spin_mutex _mutex
 
static bool _enable_print_counter = true
 

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
Date
2006

Definition at line 46 of file newton_solver.h.

Member Typedef Documentation

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

Data structure to log the information. The log is identified by the class name.

Definition at line 119 of file reference_counter.h.

Definition at line 61 of file newton_solver.h.

The type of system

Definition at line 78 of file diff_solver.h.

Member Enumeration Documentation

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)

DIVERGED_LINEAR_SOLVER_FAILURE 

The linear solver used by the DiffSolver failed to find a solution.

Definition at line 223 of file diff_solver.h.

Constructor & Destructor Documentation

libMesh::NewtonSolver::NewtonSolver ( sys_type system)
explicit

Constructor. Requires a reference to the system to be solved.

Definition at line 237 of file newton_solver.C.

238  : Parent(s),
241  brent_line_search(true),
243  minsteplength(1e-5),
246 {
247 }
std::unique_ptr< LinearSolver< Number > > _linear_solver
static std::unique_ptr< LinearSolver< Number > > build(const libMesh::Parallel::Communicator &comm_in, const SolverPackage solver_package=libMesh::default_solver_package())
libMesh::NewtonSolver::~NewtonSolver ( )
virtual

Destructor.

Definition at line 251 of file newton_solver.C.

252 {
253 }

Member Function Documentation

std::unique_ptr< DiffSolver > libMesh::DiffSolver::build ( sys_type s)
staticinherited

Factory method. Requires a reference to the system to be solved.

Returns
A NewtonSolver by default.

Definition at line 53 of file diff_solver.C.

Referenced by libMesh::TimeSolver::init().

54 {
55  return libmesh_make_unique<NewtonSolver>(s);
56 }
const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
A reference to the Parallel::Communicator object used by this mesh.

Definition at line 87 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_fd_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_mffd_residual(), libMesh::__libmesh_petsc_snes_postcheck(), libMesh::__libmesh_petsc_snes_residual(), libMesh::__libmesh_tao_equality_constraints(), libMesh::__libmesh_tao_equality_constraints_jacobian(), libMesh::__libmesh_tao_gradient(), libMesh::__libmesh_tao_hessian(), libMesh::__libmesh_tao_inequality_constraints(), libMesh::__libmesh_tao_inequality_constraints_jacobian(), libMesh::__libmesh_tao_objective(), libMesh::MeshRefinement::_coarsen_elements(), libMesh::ExactSolution::_compute_error(), libMesh::ParmetisPartitioner::_do_repartition(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add(), libMesh::EquationSystems::_read_impl(), libMesh::MeshRefinement::_refine_elements(), libMesh::MeshRefinement::_smooth_flags(), libMesh::ImplicitSystem::add_matrix(), libMesh::System::add_vector(), libMesh::EigenSparseLinearSolver< T >::adjoint_solve(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshTools::Modification::all_tri(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assemble_qoi(), libMesh::MeshCommunication::assign_global_indices(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::DofMap::attach_matrix(), libMesh::Parallel::BinSorter< KeyType, IdxType >::binsort(), libMesh::Parallel::Sort< KeyType, IdxType >::binsort(), libMesh::MeshCommunication::broadcast(), libMesh::SparseMatrix< T >::build(), libMesh::MeshTools::Generation::build_extrusion(), libMesh::Parallel::Histogram< KeyType, IdxType >::build_histogram(), libMesh::PetscNonlinearSolver< T >::build_mat_null_space(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::MeshBase::cache_elem_dims(), libMesh::System::calculate_norm(), libMesh::DofMap::check_dirichlet_bcid_consistency(), libMesh::DistributedVector< T >::clone(), libMesh::EigenSparseVector< T >::clone(), libMesh::LaspackVector< T >::clone(), libMesh::EpetraVector< T >::clone(), libMesh::PetscVector< T >::clone(), libMesh::EpetraVector< T >::close(), libMesh::Parallel::Sort< KeyType, IdxType >::communicate_bins(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshTools::create_bounding_box(), libMesh::MeshTools::create_nodal_bounding_box(), libMesh::MeshRefinement::create_parent_error_vector(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::MeshTools::create_subdomain_bounding_box(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMlibMeshSetSystem_libMesh(), DMVariableBounds_libMesh(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::EpetraVector< T >::EpetraVector(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::MeshCommunication::gather(), libMesh::MeshCommunication::gather_neighboring_elements(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::DofMap::get_info(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::EquationSystems::get_solution(), libMesh::LocationMap< T >::init(), libMesh::TopologyMap::init(), libMesh::TimeSolver::init(), libMesh::TaoOptimizationSolver< T >::init(), libMesh::PetscNonlinearSolver< T >::init(), libMesh::DistributedVector< T >::init(), libMesh::EpetraVector< T >::init(), libMesh::PetscVector< T >::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::EigenSystem::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ParmetisPartitioner::initialize(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_flags(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_p_levels(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshTools::libmesh_assert_valid_unique_ids(), libMesh::libmesh_petsc_snes_residual_helper(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::limit_overrefined_boundary(), libMesh::MeshRefinement::limit_underrefined_boundary(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::MeshCommunication::make_new_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_new_nodes_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_unique_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshCommunication::make_p_levels_parallel_consistent(), libMesh::MeshRefinement::make_refinement_compatible(), libMesh::DistributedVector< T >::max(), libMesh::FEMSystem::mesh_position_set(), libMesh::MeshSerializer::MeshSerializer(), libMesh::DistributedVector< T >::min(), libMesh::DistributedMesh::n_active_elem(), libMesh::MeshTools::n_active_levels(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::CondensedEigenSystem::n_global_non_condensed_dofs(), libMesh::MeshTools::n_levels(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::MeshTools::n_p_levels(), libMesh::BoundaryInfo::n_shellface_conds(), libMesh::DistributedMesh::parallel_max_elem_id(), libMesh::DistributedMesh::parallel_max_node_id(), libMesh::ReplicatedMesh::parallel_max_unique_id(), libMesh::DistributedMesh::parallel_max_unique_id(), libMesh::DistributedMesh::parallel_n_elem(), libMesh::DistributedMesh::parallel_n_nodes(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::MeshTools::paranoid_n_levels(), libMesh::Partitioner::partition(), libMesh::LinearPartitioner::partition_range(), libMesh::MetisPartitioner::partition_range(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::petsc_auto_fieldsplit(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshBase::prepare_for_use(), libMesh::SparseMatrix< T >::print(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::MeshCommunication::redistribute(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::CheckpointIO::select_split_config(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Partitioner::set_parent_processor_ids(), libMesh::PetscDiffSolver::setup_petsc_data(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), libMesh::split_mesh(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::Parallel::sync_element_data_by_parent_id(), libMesh::Parallel::sync_node_data_by_element_id(), libMesh::MeshRefinement::test_level_one(), libMesh::MeshRefinement::test_unflagged(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::DistributedVector< T >::zero_clone(), libMesh::EigenSparseVector< T >::zero_clone(), libMesh::LaspackVector< T >::zero_clone(), libMesh::EpetraVector< T >::zero_clone(), and libMesh::PetscVector< T >::zero_clone().

88  { return _communicator; }
const Parallel::Communicator & _communicator
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

Methods to enable/disable the reference counter output from print_info()

Definition at line 101 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::ReferenceCounter::n_objects().

102 {
103  _enable_print_counter = true;
104  return;
105 }
std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

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

Referenced by libMesh::ReferenceCounter::print_info().

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
LinearSolver<Number>& libMesh::NewtonSolver::get_linear_solver ( )
inline

Definition at line 81 of file newton_solver.h.

References _linear_solver.

82  { libmesh_assert(_linear_solver);
83  return *_linear_solver;
84  }
std::unique_ptr< LinearSolver< Number > > _linear_solver
const LinearSolver<Number>& libMesh::NewtonSolver::get_linear_solver ( ) const
inline

Definition at line 86 of file newton_solver.h.

References _linear_solver.

87  { libmesh_assert(_linear_solver);
88  return *_linear_solver;
89  }
std::unique_ptr< LinearSolver< Number > > _linear_solver
void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 185 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

186 {
187  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188  std::pair<unsigned int, unsigned int> & p = _counts[name];
189 
190  p.first++;
191 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 198 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

199 {
200  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
201  std::pair<unsigned int, unsigned int> & p = _counts[name];
202 
203  p.second++;
204 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
void libMesh::NewtonSolver::init ( )
virtual

The initialization function. This method is used to initialize internal data structures before a simulation begins.

Reimplemented from libMesh::DiffSolver.

Definition at line 257 of file newton_solver.C.

References _linear_solver, libMesh::DiffSolver::_system, libMesh::DiffSolver::init(), libMesh::System::name(), and libMesh::on_command_line().

258 {
259  Parent::init();
260 
261  if (libMesh::on_command_line("--solver_system_names"))
262  _linear_solver->init((_system.name()+"_").c_str());
263  else
264  _linear_solver->init();
265 
266  _linear_solver->init_names(_system);
267 }
std::unique_ptr< LinearSolver< Number > > _linear_solver
const std::string & name() const
Definition: system.h:1998
virtual void init()
Definition: diff_solver.C:69
sys_type & _system
Definition: diff_solver.h:319
bool on_command_line(const std::string &arg)
Definition: libmesh.C:920
Real libMesh::NewtonSolver::line_search ( Real  tol,
Real  last_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.

Returns
The substep size.

Definition at line 38 of file newton_solver.C.

References libMesh::DiffSolver::_outer_iterations, libMesh::DiffSolver::_solve_result, libMesh::DiffSolver::_system, std::abs(), libMesh::NumericVector< T >::add(), libMesh::ImplicitSystem::assembly(), brent_line_search, libMesh::NumericVector< T >::close(), libMesh::DiffSolver::continue_after_backtrack_failure, libMesh::DiffSolver::DIVERGED_BACKTRACKING_FAILURE, libMesh::NumericVector< T >::l2_norm(), libMesh::libmesh_isnan(), std::max(), std::min(), minsteplength, libMesh::out, libMesh::DiffSolver::quiet, libMesh::Real, require_finite_residual, require_residual_reduction, libMesh::ExplicitSystem::rhs, libMesh::SIGN(), libMesh::System::update(), and libMesh::DiffSolver::verbose.

Referenced by solve().

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
52  NumericVector<Number> & rhs = *(_system.rhs);
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  // We may need to localize a parallel solution
83  _system.update();
84 
85  // Check residual with fractional Newton step
86  _system.assembly (true, false);
87 
88  rhs.close();
89  current_residual = rhs.l2_norm();
90  if (verbose)
91  libMesh::out << " Current Residual: "
92  << current_residual << std::endl;
93 
94  if (bx/2. < minsteplength &&
95  (libmesh_isnan(current_residual) ||
96  (current_residual > last_residual)))
97  {
98  libMesh::out << "Inexact Newton step FAILED at step "
99  << _outer_iterations << std::endl;
100 
102  {
103  libmesh_convergence_failure();
104  }
105  else
106  {
107  libMesh::out << "Continuing anyway ..." << std::endl;
109  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  if (!brent_line_search)
118  return bx;
119 
120  // Brent's method adapted from Numerical Recipes in C, ch. 10.2
121  Real e = 0.;
122 
123  Real x = bx, w = bx, v = bx;
124 
125  // Residuals at bx
126  Real fx = current_residual,
127  fw = current_residual,
128  fv = current_residual;
129 
130  // Max iterations for Brent's method loop
131  const unsigned int max_i = 20;
132 
133  // for golden ratio steps
134  const Real golden_ratio = 1.-(std::sqrt(5.)-1.)/2.;
135 
136  for (unsigned int i=1; i <= max_i; i++)
137  {
138  Real xm = (ax+cx)*0.5;
139  Real tol1 = tol * std::abs(x) + tol*tol;
140  Real tol2 = 2.0 * tol1;
141 
142  // Test if we're done
143  if (std::abs(x-xm) <= (tol2 - 0.5 * (cx - ax)))
144  return x;
145 
146  Real d;
147 
148  // Construct a parabolic fit
149  if (std::abs(e) > tol1)
150  {
151  Real r = (x-w)*(fx-fv);
152  Real q = (x-v)*(fx-fw);
153  Real p = (x-v)*q-(x-w)*r;
154  q = 2. * (q-r);
155  if (q > 0.)
156  p = -p;
157  else
158  q = std::abs(q);
159  if (std::abs(p) >= std::abs(0.5*q*e) ||
160  p <= q * (ax-x) ||
161  p >= q * (cx-x))
162  {
163  // Take a golden section step
164  e = x >= xm ? ax-x : cx-x;
165  d = golden_ratio * e;
166  }
167  else
168  {
169  // Take a parabolic fit step
170  d = p/q;
171  if (x+d-ax < tol2 || cx-(x+d) < tol2)
172  d = SIGN(tol1, xm - x);
173  }
174  }
175  else
176  {
177  // Take a golden section step
178  e = x >= xm ? ax-x : cx-x;
179  d = golden_ratio * e;
180  }
181 
182  Real u = std::abs(d) >= tol1 ? x+d : x + SIGN(tol1,d);
183 
184  // Assemble the residual at the new steplength u
185  newton_iterate.add (bx - u, linear_solution);
186  newton_iterate.close();
187  bx = u;
188  if (verbose)
189  libMesh::out << " Shrinking Newton step to "
190  << bx << std::endl;
191 
192  // We may need to localize a parallel solution
193  _system.update();
194  _system.assembly (true, false);
195 
196  rhs.close();
197  Real fu = current_residual = rhs.l2_norm();
198  if (verbose)
199  libMesh::out << " Current Residual: "
200  << fu << std::endl;
201 
202  if (fu <= fx)
203  {
204  if (u >= x)
205  ax = x;
206  else
207  cx = x;
208  v = w; w = x; x = u;
209  fv = fw; fw = fx; fx = fu;
210  }
211  else
212  {
213  if (u < x)
214  ax = u;
215  else
216  cx = u;
217  if (fu <= fw || w == x)
218  {
219  v = w; w = u;
220  fv = fw; fw = fu;
221  }
222  else if (fu <= fv || v == x || v == w)
223  {
224  v = u;
225  fv = fu;
226  }
227  }
228  }
229 
230  if (!quiet)
231  libMesh::out << "Warning! Too many iterations used in Brent line search!"
232  << std::endl;
233  return bx;
234 }
double abs(double a)
NumericVector< Number > * rhs
long double max(long double a, double b)
unsigned int _solve_result
Definition: diff_solver.h:327
T SIGN(T a, T b)
Definition: newton_solver.C:33
virtual void assembly(bool, bool, bool=false, bool=false)
sys_type & _system
Definition: diff_solver.h:319
virtual void close()=0
virtual void update()
Definition: system.C:406
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool continue_after_backtrack_failure
Definition: diff_solver.h:181
bool libmesh_isnan(float a)
unsigned int _outer_iterations
Definition: diff_solver.h:309
virtual void add(const numeric_index_type i, const T value)=0
OStreamProxy out(std::cout)
long double min(long double a, double b)
static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited
processor_id_type libMesh::ParallelObject::n_processors ( ) const
inlineinherited
Returns
The number of processors in the group.

Definition at line 93 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::size().

Referenced by libMesh::ParmetisPartitioner::_do_repartition(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::DistributedMesh::add_elem(), libMesh::DistributedMesh::add_node(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assembly(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::AztecLinearSolver< T >::AztecLinearSolver(), libMesh::MeshCommunication::broadcast(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::DistributedMesh::clear(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::UnstructuredMesh::create_pid_mesh(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::DistributedMesh::DistributedMesh(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshCommunication::gather(), libMesh::MeshCommunication::gather_neighboring_elements(), libMesh::MeshBase::get_info(), libMesh::EquationSystems::get_solution(), libMesh::DistributedVector< T >::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::DistributedMesh::insert_elem(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::DofMap::local_variable_indices(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshBase::n_active_elem_on_proc(), libMesh::MeshBase::n_elem_on_proc(), libMesh::MeshBase::n_nodes_on_proc(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::Partitioner::partition(), libMesh::MeshBase::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::PetscLinearSolver< T >::PetscLinearSolver(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::SparseMatrix< T >::print(), libMesh::NameBasedIO::read(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::CheckpointIO::read_connectivity(), libMesh::XdrIO::read_header(), libMesh::CheckpointIO::read_nodes(), libMesh::MeshCommunication::redistribute(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::Partitioner::repartition(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

94  { return cast_int<processor_id_type>(_communicator.size()); }
unsigned int size() const
Definition: parallel.h:725
const Parallel::Communicator & _communicator
void libMesh::NewtonSolver::print_convergence ( unsigned int  step_num,
Real  current_residual,
Real  step_norm,
bool  linear_solve_finished 
)
protected

This prints output for the convergence criteria based on by the given residual and step size.

Definition at line 618 of file newton_solver.C.

References libMesh::DiffSolver::absolute_residual_tolerance, libMesh::DiffSolver::absolute_step_tolerance, libMesh::DiffSolver::max_residual_norm, libMesh::DiffSolver::max_solution_norm, libMesh::out, libMesh::DiffSolver::relative_residual_tolerance, libMesh::DiffSolver::relative_step_tolerance, and libMesh::DiffSolver::verbose.

Referenced by solve().

622 {
623  // Is our absolute residual low enough?
624  if (current_residual < absolute_residual_tolerance)
625  {
626  libMesh::out << " Nonlinear solver converged, step " << step_num
627  << ", residual " << current_residual
628  << std::endl;
629  }
631  {
632  if (verbose)
633  libMesh::out << " Nonlinear solver current_residual "
634  << current_residual << " > "
635  << (absolute_residual_tolerance) << std::endl;
636  }
637 
638  // Is our relative residual low enough?
639  if ((current_residual / max_residual_norm) <
641  {
642  libMesh::out << " Nonlinear solver converged, step " << step_num
643  << ", residual reduction "
644  << current_residual / max_residual_norm
645  << " < " << relative_residual_tolerance
646  << std::endl;
647  }
649  {
650  if (verbose)
651  libMesh::out << " Nonlinear solver relative residual "
652  << (current_residual / max_residual_norm)
653  << " > " << relative_residual_tolerance
654  << std::endl;
655  }
656 
657  // For incomplete linear solves, it's not safe to test step sizes
658  if (!linear_solve_finished)
659  return;
660 
661  // Is our absolute Newton step size small enough?
662  if (step_norm < absolute_step_tolerance)
663  {
664  libMesh::out << " Nonlinear solver converged, step " << step_num
665  << ", absolute step size "
666  << step_norm
667  << " < " << absolute_step_tolerance
668  << std::endl;
669  }
670  else if (absolute_step_tolerance)
671  {
672  if (verbose)
673  libMesh::out << " Nonlinear solver absolute step size "
674  << step_norm
675  << " > " << absolute_step_tolerance
676  << std::endl;
677  }
678 
679  // Is our relative Newton step size small enough?
680  if (step_norm / max_solution_norm <
682  {
683  libMesh::out << " Nonlinear solver converged, step " << step_num
684  << ", relative step size "
685  << (step_norm / max_solution_norm)
686  << " < " << relative_step_tolerance
687  << std::endl;
688  }
689  else if (relative_step_tolerance)
690  {
691  if (verbose)
692  libMesh::out << " Nonlinear solver relative step size "
693  << (step_norm / max_solution_norm)
694  << " > " << relative_step_tolerance
695  << std::endl;
696  }
697 }
Real absolute_residual_tolerance
Definition: diff_solver.h:192
OStreamProxy out(std::cout)
Real relative_residual_tolerance
Definition: diff_solver.h:193
void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

Referenced by libMesh::LibMeshInit::LibMeshInit().

89 {
91  out_stream << ReferenceCounter::get_info();
92 }
static std::string get_info()
processor_id_type libMesh::ParallelObject::processor_id ( ) const
inlineinherited
Returns
The rank of this processor in the group.

Definition at line 99 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::rank().

Referenced by libMesh::BoundaryInfo::_find_id_maps(), libMesh::EquationSystems::_read_impl(), libMesh::DistributedMesh::add_elem(), libMesh::BoundaryInfo::add_elements(), libMesh::DofMap::add_neighbors_to_send_list(), libMesh::DistributedMesh::add_node(), libMesh::MeshRefinement::add_node(), libMesh::MeshTools::Modification::all_tri(), libMesh::FEMSystem::assembly(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::MeshCommunication::broadcast(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::Nemesis_IO_Helper::build_element_and_node_maps(), libMesh::ParmetisPartitioner::build_graph(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::DofMap::build_sparsity(), libMesh::DistributedMesh::clear(), libMesh::ExodusII_IO_Helper::close(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::compute_communication_map_parameters(), libMesh::Nemesis_IO_Helper::compute_internal_and_border_elems_and_internal_nodes(), libMesh::Nemesis_IO_Helper::compute_node_communication_maps(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::ExodusII_IO_Helper::create(), libMesh::DistributedMesh::delete_elem(), libMesh::DistributedMesh::delete_node(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::DistributedMesh::DistributedMesh(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshFunction::find_element(), libMesh::MeshFunction::find_elements(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::MeshCommunication::gather(), libMesh::MeshCommunication::gather_neighboring_elements(), libMesh::Nemesis_IO_Helper::get_cmap_params(), libMesh::Nemesis_IO_Helper::get_eb_info_global(), libMesh::Nemesis_IO_Helper::get_elem_cmap(), libMesh::Nemesis_IO_Helper::get_elem_map(), libMesh::MeshBase::get_info(), libMesh::DofMap::get_info(), libMesh::Nemesis_IO_Helper::get_init_global(), libMesh::Nemesis_IO_Helper::get_init_info(), libMesh::Nemesis_IO_Helper::get_loadbal_param(), libMesh::Nemesis_IO_Helper::get_node_cmap(), libMesh::Nemesis_IO_Helper::get_node_map(), libMesh::Nemesis_IO_Helper::get_ns_param_global(), libMesh::EquationSystems::get_solution(), libMesh::Nemesis_IO_Helper::get_ss_param_global(), libMesh::SparsityPattern::Build::handle_vi_vj(), libMesh::DistributedVector< T >::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::ExodusII_IO_Helper::initialize(), libMesh::ExodusII_IO_Helper::initialize_element_variables(), libMesh::ExodusII_IO_Helper::initialize_global_variables(), libMesh::ExodusII_IO_Helper::initialize_nodal_variables(), libMesh::DistributedMesh::insert_elem(), libMesh::DofMap::is_evaluable(), libMesh::SparsityPattern::Build::join(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), libMesh::DofMap::local_variable_indices(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshBase::n_active_local_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::System::n_local_dofs(), libMesh::MeshBase::n_local_elem(), libMesh::MeshBase::n_local_nodes(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::BoundaryInfo::n_shellface_conds(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SparsityPattern::Build::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::MetisPartitioner::partition_range(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::SparseMatrix< T >::print(), libMesh::NumericVector< T >::print_global(), libMesh::Nemesis_IO_Helper::put_cmap_params(), libMesh::Nemesis_IO_Helper::put_elem_cmap(), libMesh::Nemesis_IO_Helper::put_elem_map(), libMesh::Nemesis_IO_Helper::put_loadbal_param(), libMesh::Nemesis_IO_Helper::put_node_cmap(), libMesh::Nemesis_IO_Helper::put_node_map(), libMesh::NameBasedIO::read(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read(), libMesh::ExodusII_IO_Helper::read_elem_num_map(), libMesh::ExodusII_IO_Helper::read_global_values(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::MeshCommunication::redistribute(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::CheckpointIO::select_split_config(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::Parallel::Packing< Node * >::unpack(), libMesh::Parallel::Packing< Elem * >::unpack(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::MeshTools::weight(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::CheckpointIO::write(), libMesh::EquationSystems::write(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO_Helper::write_element_values(), libMesh::ExodusII_IO_Helper::write_elements(), libMesh::ExodusII_IO::write_global_data(), libMesh::ExodusII_IO_Helper::write_global_values(), libMesh::ExodusII_IO::write_information_records(), libMesh::ExodusII_IO_Helper::write_information_records(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates(), libMesh::UCDIO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::ExodusII_IO_Helper::write_nodal_values(), libMesh::Nemesis_IO_Helper::write_nodesets(), libMesh::ExodusII_IO_Helper::write_nodesets(), libMesh::XdrIO::write_serialized_bc_names(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::Nemesis_IO_Helper::write_sidesets(), libMesh::ExodusII_IO_Helper::write_sidesets(), libMesh::ExodusII_IO::write_timestep(), libMesh::ExodusII_IO_Helper::write_timestep(), and libMesh::ExodusII_IO::write_timestep_discontinuous().

100  { return cast_int<processor_id_type>(_communicator.rank()); }
const Parallel::Communicator & _communicator
unsigned int rank() const
Definition: parallel.h:723
void libMesh::NewtonSolver::reinit ( )
virtual

The reinitialization function. This method is used after changes in the mesh.

Reimplemented from libMesh::DiffSolver.

Definition at line 271 of file newton_solver.C.

References _linear_solver, libMesh::DiffSolver::_system, and libMesh::DiffSolver::reinit().

272 {
273  Parent::reinit();
274 
275  _linear_solver->clear();
276 
277  _linear_solver->init_names(_system);
278 }
std::unique_ptr< LinearSolver< Number > > _linear_solver
virtual void reinit()
Definition: diff_solver.C:60
sys_type & _system
Definition: diff_solver.h:319
unsigned int libMesh::NewtonSolver::solve ( )
virtual

This method performs a solve, using an inexact Newton-Krylov method with line search.

Implements libMesh::DiffSolver.

Definition at line 282 of file newton_solver.C.

References libMesh::DiffSolver::_inner_iterations, _linear_solver, libMesh::DiffSolver::_outer_iterations, libMesh::DiffSolver::_solve_result, libMesh::DiffSolver::_system, libMesh::DiffSolver::absolute_residual_tolerance, libMesh::DiffSolver::absolute_step_tolerance, libMesh::NumericVector< T >::add(), libMesh::ImplicitSystem::assembly(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::DiffSolver::continue_after_max_iterations, libMesh::DiffSolver::CONVERGED_ABSOLUTE_RESIDUAL, libMesh::DiffSolver::CONVERGED_ABSOLUTE_STEP, libMesh::DiffSolver::CONVERGED_RELATIVE_RESIDUAL, libMesh::DiffSolver::CONVERGED_RELATIVE_STEP, libMesh::DiffSolver::DIVERGED_BACKTRACKING_FAILURE, libMesh::DiffSolver::DIVERGED_LINEAR_SOLVER_FAILURE, libMesh::DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS, libMesh::DofMap::enforce_constraints_exactly(), libMesh::System::get_dof_map(), libMesh::DiffSolver::initial_linear_tolerance, libMesh::DiffSolver::INVALID_SOLVE_RESULT, libMesh::NumericVector< T >::l2_norm(), libMesh::libmesh_isnan(), line_search(), libMesh::DiffSolver::linear_solution_monitor, linear_tolerance_multiplier, libMesh::ImplicitSystem::matrix, std::max(), libMesh::DiffSolver::max_linear_iterations, libMesh::DiffSolver::max_nonlinear_iterations, libMesh::DiffSolver::max_residual_norm, libMesh::DiffSolver::max_solution_norm, std::min(), libMesh::DiffSolver::minimum_linear_tolerance, libMesh::out, print_convergence(), libMesh::DiffSolver::quiet, libMesh::Real, libMesh::DiffSolver::relative_residual_tolerance, libMesh::DiffSolver::relative_step_tolerance, libMesh::ImplicitSystem::request_matrix(), require_finite_residual, require_residual_reduction, libMesh::ExplicitSystem::rhs, libMesh::System::solution, test_convergence(), libMesh::TOLERANCE, track_linear_convergence, libMesh::System::update(), libMesh::DiffSolver::verbose, libMesh::NumericVector< T >::zero(), and libMesh::NumericVector< T >::zero_clone().

283 {
284  LOG_SCOPE("solve()", "NewtonSolver");
285 
286  // Reset any prior solve result
288 
289  NumericVector<Number> & newton_iterate = *(_system.solution);
290 
291  std::unique_ptr<NumericVector<Number>> linear_solution_ptr = newton_iterate.zero_clone();
292  NumericVector<Number> & linear_solution = *linear_solution_ptr;
293  NumericVector<Number> & rhs = *(_system.rhs);
294 
295  newton_iterate.close();
296  linear_solution.close();
297  rhs.close();
298 
299 #ifdef LIBMESH_ENABLE_CONSTRAINTS
301 #endif
302 
303  SparseMatrix<Number> & matrix = *(_system.matrix);
304 
305  // Set starting linear tolerance
306  Real current_linear_tolerance = initial_linear_tolerance;
307 
308  // Start counting our linear solver steps
309  _inner_iterations = 0;
310 
311  // Now we begin the nonlinear loop
314  {
315  // We may need to localize a parallel solution
316  _system.update();
317 
318  if (verbose)
319  libMesh::out << "Assembling the System" << std::endl;
320 
321  _system.assembly(true, true);
322  rhs.close();
323  Real current_residual = rhs.l2_norm();
324 
325  if (libmesh_isnan(current_residual))
326  {
327  libMesh::out << " Nonlinear solver DIVERGED at step "
329  << " with residual Not-a-Number"
330  << std::endl;
331  libmesh_convergence_failure();
332  continue;
333  }
334 
335  if (current_residual <= absolute_residual_tolerance)
336  {
337  if (verbose)
338  libMesh::out << "Linear solve unnecessary; residual "
339  << current_residual
340  << " meets absolute tolerance "
342  << std::endl;
343 
344  // We're not doing a solve, but other code may reuse this
345  // matrix.
346  matrix.close();
347 
349  if (current_residual == 0)
350  {
353  if (absolute_step_tolerance > 0)
355  if (relative_step_tolerance > 0)
357  }
358 
359  break;
360  }
361 
362  // Prepare to take incomplete steps
363  Real last_residual = current_residual;
364 
365  max_residual_norm = std::max (current_residual,
367 
368  // Compute the l2 norm of the whole solution
369  Real norm_total = newton_iterate.l2_norm();
370 
372 
373  if (verbose)
374  libMesh::out << "Nonlinear Residual: "
375  << current_residual << std::endl;
376 
377  // Make sure our linear tolerance is low enough
378  current_linear_tolerance = std::min (current_linear_tolerance,
379  current_residual * linear_tolerance_multiplier);
380 
381  // But don't let it be too small
382  if (current_linear_tolerance < minimum_linear_tolerance)
383  {
384  current_linear_tolerance = minimum_linear_tolerance;
385  }
386 
387  // If starting the nonlinear solve with a really good initial guess, we dont want to set an absurd linear tolerance
388  current_linear_tolerance = std::max(current_linear_tolerance, absolute_residual_tolerance / current_residual / 10.0);
389 
390  // At this point newton_iterate is the current guess, and
391  // linear_solution is now about to become the NEGATIVE of the next
392  // Newton step.
393 
394  // Our best initial guess for the linear_solution is zero!
395  linear_solution.zero();
396 
397  if (verbose)
398  libMesh::out << "Linear solve starting, tolerance "
399  << current_linear_tolerance << std::endl;
400 
401  // Solve the linear system.
402  const std::pair<unsigned int, Real> rval =
403  _linear_solver->solve (matrix, _system.request_matrix("Preconditioner"),
404  linear_solution, rhs, current_linear_tolerance,
406 
408  {
409  LinearConvergenceReason linear_c_reason = _linear_solver->get_converged_reason();
410 
411  // Check if something went wrong during the linear solve
412  if (linear_c_reason < 0)
413  {
414  // The linear solver failed somehow
416  // Print a message
417  libMesh::out << "Linear solver failed during Newton step, dropping out."
418  << std::endl;
419  break;
420  }
421  }
422 
423  // We may need to localize a parallel solution
424  _system.update ();
425  // The linear solver may not have fit our constraints exactly
426 #ifdef LIBMESH_ENABLE_CONSTRAINTS
428  (_system, &linear_solution, /* homogeneous = */ true);
429 #endif
430 
431  const unsigned int linear_steps = rval.first;
432  libmesh_assert_less_equal (linear_steps, max_linear_iterations);
433  _inner_iterations += linear_steps;
434 
435  const bool linear_solve_finished =
436  !(linear_steps == max_linear_iterations);
437 
438  if (verbose)
439  libMesh::out << "Linear solve finished, step " << linear_steps
440  << ", residual " << rval.second
441  << std::endl;
442 
443  // Compute the l2 norm of the nonlinear update
444  Real norm_delta = linear_solution.l2_norm();
445 
446  if (verbose)
447  libMesh::out << "Trying full Newton step" << std::endl;
448  // Take a full Newton step
449  newton_iterate.add (-1., linear_solution);
450  newton_iterate.close();
451 
452  if (this->linear_solution_monitor.get())
453  {
454  // Compute the l2 norm of the whole solution
455  norm_total = newton_iterate.l2_norm();
456  rhs.close();
457  (*this->linear_solution_monitor)(linear_solution, norm_delta,
458  newton_iterate, norm_total,
459  rhs, rhs.l2_norm(), _outer_iterations);
460  }
461 
462  // Check residual with full Newton step, if that's useful for determining
463  // whether to line search, whether to quit early, or whether to die after
464  // hitting our max iteration count
465  if (this->require_residual_reduction ||
466  this->require_finite_residual ||
467  _outer_iterations+1 < max_nonlinear_iterations ||
469  {
470  _system.update ();
471  _system.assembly(true, false);
472 
473  rhs.close();
474  current_residual = rhs.l2_norm();
475  if (verbose)
476  libMesh::out << " Current Residual: "
477  << current_residual << std::endl;
478 
479  // don't fiddle around if we've already converged
480  if (test_convergence(current_residual, norm_delta,
481  linear_solve_finished &&
482  current_residual <= last_residual))
483  {
484  if (!quiet)
485  print_convergence(_outer_iterations, current_residual,
486  norm_delta, linear_solve_finished &&
487  current_residual <= last_residual);
489  break; // out of _outer_iterations for loop
490  }
491  }
492 
493  // since we're not converged, backtrack if necessary
494  Real steplength =
495  this->line_search(std::sqrt(TOLERANCE),
496  last_residual, current_residual,
497  newton_iterate, linear_solution);
498  norm_delta *= steplength;
499 
500  // Check to see if backtracking failed,
501  // and break out of the nonlinear loop if so...
503  {
505  break; // out of _outer_iterations for loop
506  }
507 
508  if (_outer_iterations + 1 >= max_nonlinear_iterations)
509  {
510  libMesh::out << " Nonlinear solver reached maximum step "
511  << max_nonlinear_iterations << ", latest evaluated residual "
512  << current_residual << std::endl;
514  {
516  libMesh::out << " Continuing..." << std::endl;
517  }
518  else
519  {
520  libmesh_convergence_failure();
521  }
522  continue;
523  }
524 
525  // Compute the l2 norm of the whole solution
526  norm_total = newton_iterate.l2_norm();
527 
529 
530  // Print out information for the
531  // nonlinear iterations.
532  if (verbose)
533  libMesh::out << " Nonlinear step: |du|/|u| = "
534  << norm_delta / norm_total
535  << ", |du| = " << norm_delta
536  << std::endl;
537 
538  // Terminate the solution iteration if the difference between
539  // this iteration and the last is sufficiently small.
540  if (test_convergence(current_residual, norm_delta / steplength,
541  linear_solve_finished))
542  {
543  if (!quiet)
544  print_convergence(_outer_iterations, current_residual,
545  norm_delta / steplength,
546  linear_solve_finished);
548  break; // out of _outer_iterations for loop
549  }
550  } // end nonlinear loop
551 
552  // The linear solver may not have fit our constraints exactly
553 #ifdef LIBMESH_ENABLE_CONSTRAINTS
555 #endif
556 
557  // We may need to localize a parallel solution
558  _system.update ();
559 
560  // Make sure we are returning something sensible as the
561  // _solve_result, except in the edge case where we weren't really asked to
562  // solve.
563  libmesh_assert (_solve_result != DiffSolver::INVALID_SOLVE_RESULT ||
564  !max_nonlinear_iterations);
565 
566  return _solve_result;
567 }
bool continue_after_max_iterations
Definition: diff_solver.h:175
std::unique_ptr< LinearSolver< Number > > _linear_solver
unsigned int max_nonlinear_iterations
Definition: diff_solver.h:157
Real absolute_residual_tolerance
Definition: diff_solver.h:192
NumericVector< Number > * rhs
static const Real TOLERANCE
long double max(long double a, double b)
unsigned int _solve_result
Definition: diff_solver.h:327
unsigned int _inner_iterations
Definition: diff_solver.h:314
virtual void assembly(bool, bool, bool=false, bool=false)
const DofMap & get_dof_map() const
Definition: system.h:2030
sys_type & _system
Definition: diff_solver.h:319
unsigned int max_linear_iterations
Definition: diff_solver.h:149
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
const SparseMatrix< Number > * request_matrix(const std::string &mat_name) const
virtual void update()
Definition: system.C:406
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool libmesh_isnan(float a)
unsigned int _outer_iterations
Definition: diff_solver.h:309
SparseMatrix< Number > * matrix
bool test_convergence(Real current_residual, Real step_norm, bool linear_solve_finished)
Real line_search(Real tol, Real last_residual, Real &current_residual, NumericVector< Number > &newton_iterate, const NumericVector< Number > &linear_solution)
Definition: newton_solver.C:38
void print_convergence(unsigned int step_num, Real current_residual, Real step_norm, bool linear_solve_finished)
OStreamProxy out(std::cout)
long double min(long double a, double b)
Real relative_residual_tolerance
Definition: diff_solver.h:193
std::unique_ptr< LinearSolutionMonitor > linear_solution_monitor
Definition: diff_solver.h:289
unsigned int libMesh::DiffSolver::solve_result ( )
inlineinherited
Returns
The value of the SolveResult from the last solve.

Definition at line 133 of file diff_solver.h.

133 { return _solve_result; }
unsigned int _solve_result
Definition: diff_solver.h:327
const sys_type& libMesh::DiffSolver::system ( ) const
inlineinherited
Returns
A constant reference to the system we are solving.

Definition at line 138 of file diff_solver.h.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), and libMesh::__libmesh_petsc_diff_solver_residual().

138 { return _system; }
sys_type & _system
Definition: diff_solver.h:319
sys_type& libMesh::DiffSolver::system ( )
inlineinherited
Returns
A writable reference to the system we are solving.

Definition at line 143 of file diff_solver.h.

143 { return _system; }
sys_type & _system
Definition: diff_solver.h:319
bool libMesh::NewtonSolver::test_convergence ( Real  current_residual,
Real  step_norm,
bool  linear_solve_finished 
)
protected
Returns
true if a convergence criterion has been passed by the given residual and step size; false otherwise.

Definition at line 571 of file newton_solver.C.

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

Referenced by solve().

574 {
575  // We haven't converged unless we pass a convergence test
576  bool has_converged = false;
577 
578  // Is our absolute residual low enough?
579  if (current_residual < absolute_residual_tolerance)
580  {
582  has_converged = true;
583  }
584 
585  // Is our relative residual low enough?
586  if ((current_residual / max_residual_norm) <
588  {
590  has_converged = true;
591  }
592 
593  // For incomplete linear solves, it's not safe to test step sizes
594  if (!linear_solve_finished)
595  {
596  return has_converged;
597  }
598 
599  // Is our absolute Newton step size small enough?
600  if (step_norm < absolute_step_tolerance)
601  {
603  has_converged = true;
604  }
605 
606  // Is our relative Newton step size small enough?
607  if (step_norm / max_solution_norm <
609  {
611  has_converged = true;
612  }
613 
614  return has_converged;
615 }
Real absolute_residual_tolerance
Definition: diff_solver.h:192
unsigned int _solve_result
Definition: diff_solver.h:327
Real relative_residual_tolerance
Definition: diff_solver.h:193
unsigned int libMesh::DiffSolver::total_inner_iterations ( )
inlineinherited
Returns
The number of "inner" (e.g. Krylov) iterations required by the last solve.

Definition at line 128 of file diff_solver.h.

128 { return _inner_iterations; }
unsigned int _inner_iterations
Definition: diff_solver.h:314
unsigned int libMesh::DiffSolver::total_outer_iterations ( )
inlineinherited
Returns
The number of "outer" (e.g. quasi-Newton) iterations required by the last solve.

Definition at line 122 of file diff_solver.h.

122 { return _outer_iterations; }
unsigned int _outer_iterations
Definition: diff_solver.h:309

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited
bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

Flag to control whether reference count information is printed when print_info is called.

Definition at line 143 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

unsigned int libMesh::DiffSolver::_inner_iterations
protectedinherited

The number of inner iterations used by the last solve

Definition at line 314 of file diff_solver.h.

Referenced by solve().

std::unique_ptr<LinearSolver<Number> > libMesh::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 153 of file newton_solver.h.

Referenced by get_linear_solver(), init(), reinit(), and solve().

Threads::spin_mutex libMesh::ReferenceCounter::_mutex
staticprotectedinherited

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 137 of file reference_counter.h.

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects
staticprotectedinherited

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 132 of file reference_counter.h.

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

unsigned int libMesh::DiffSolver::_outer_iterations
protectedinherited

The number of outer iterations used by the last solve

Definition at line 309 of file diff_solver.h.

Referenced by line_search(), and solve().

unsigned int libMesh::DiffSolver::_solve_result
protectedinherited

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 327 of file diff_solver.h.

Referenced by line_search(), solve(), and test_convergence().

sys_type& libMesh::DiffSolver::_system
protectedinherited

A reference to the system we are solving.

Definition at line 319 of file diff_solver.h.

Referenced by init(), line_search(), reinit(), libMesh::PetscDiffSolver::setup_petsc_data(), solve(), and libMesh::PetscDiffSolver::solve().

Real libMesh::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 192 of file diff_solver.h.

Referenced by print_convergence(), solve(), and test_convergence().

Real libMesh::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 204 of file diff_solver.h.

Referenced by print_convergence(), solve(), and test_convergence().

bool libMesh::NewtonSolver::brent_line_search

If 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 120 of file newton_solver.h.

Referenced by line_search().

bool libMesh::DiffSolver::continue_after_backtrack_failure
inherited

Defaults to false, telling the DiffSolver to throw an error when the backtracking scheme fails to find a descent direction.

Definition at line 181 of file diff_solver.h.

Referenced by line_search().

bool libMesh::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 175 of file diff_solver.h.

Referenced by solve().

Real libMesh::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 211 of file diff_solver.h.

Referenced by solve().

std::unique_ptr<LinearSolutionMonitor> libMesh::DiffSolver::linear_solution_monitor
inherited

Pointer to functor which is called right after each linear solve

Definition at line 289 of file diff_solver.h.

Referenced by libMesh::__libmesh_petsc_diff_solver_monitor(), and solve().

Real libMesh::NewtonSolver::linear_tolerance_multiplier

The 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 143 of file newton_solver.h.

Referenced by solve().

unsigned int libMesh::DiffSolver::max_linear_iterations
inherited

Each linear solver step should exit after max_linear_iterations is exceeded.

Definition at line 149 of file diff_solver.h.

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

unsigned int libMesh::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 157 of file diff_solver.h.

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

Real libMesh::DiffSolver::max_residual_norm
protectedinherited

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 304 of file diff_solver.h.

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

Real libMesh::DiffSolver::max_solution_norm
protectedinherited

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 297 of file diff_solver.h.

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

Real libMesh::DiffSolver::minimum_linear_tolerance
inherited

The tolerance for linear solves is kept above this minimum

Definition at line 216 of file diff_solver.h.

Referenced by solve().

Real libMesh::NewtonSolver::minsteplength

If the quasi-Newton step length must be reduced to below this factor to give a residual reduction, then the Newton solver dies with an error message. It is currently set to 1e-5 by default.

Definition at line 137 of file newton_solver.h.

Referenced by line_search().

bool libMesh::DiffSolver::quiet
inherited

The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is true.

Definition at line 163 of file diff_solver.h.

Referenced by line_search(), and solve().

Real libMesh::DiffSolver::relative_residual_tolerance
inherited

Definition at line 193 of file diff_solver.h.

Referenced by print_convergence(), solve(), and test_convergence().

Real libMesh::DiffSolver::relative_step_tolerance
inherited

Definition at line 205 of file diff_solver.h.

Referenced by print_convergence(), solve(), and test_convergence().

bool libMesh::NewtonSolver::require_finite_residual

If 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 an infinite or NaN residual. It is currently set to true by default; set it to false to avoid unnecessary residual assembly on well-behaved systems.

Definition at line 107 of file newton_solver.h.

Referenced by line_search(), and solve().

bool libMesh::NewtonSolver::require_residual_reduction

If 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 98 of file newton_solver.h.

Referenced by line_search(), and solve().

bool libMesh::NewtonSolver::track_linear_convergence

If set to true, check for convergence of the linear solve. If no convergence is acquired during the linear solve, the nonlinear solve fails with DiffSolver::DIVERGED_LINEAR_SOLVER_FAILURE. Enabled by default as nonlinear convergence is very difficult, if the linear solver is not converged.

Definition at line 129 of file newton_solver.h.

Referenced by solve().

bool libMesh::DiffSolver::verbose
inherited

The documentation for this class was generated from the following files: