#include <petsc_nonlinear_solver.h>

Inheritance diagram for libMesh::PetscNonlinearSolver< T >:

Public Types
typedef NonlinearImplicitSystem	sys_type
Public Member Functions
	PetscNonlinearSolver (sys_type &system)
	~PetscNonlinearSolver ()
virtual void	clear ()
virtual void	init (const char *name=NULL)
SNES	snes ()
virtual std::pair< unsigned int, Real >	solve (SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)
virtual void	print_converged_reason ()
SNESConvergedReason	get_converged_reason ()
virtual int	get_total_linear_iterations ()
virtual unsigned	get_current_nonlinear_iteration_number () const
void	set_residual_zero_out (bool state)
void	set_jacobian_zero_out (bool state)
void	use_default_monitor (bool state)
bool	initialized () const
const sys_type &	system () const
sys_type &	system ()
void	attach_preconditioner (Preconditioner< T > *preconditioner)
const Parallel::Communicator &	comm () const
processor_id_type	n_processors () const
processor_id_type	processor_id () const
Static Public Member Functions
static UniquePtr < NonlinearSolver< T > >	build (sys_type &s, const SolverPackage solver_package=libMesh::default_solver_package())
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
void(*	residual )(const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)
NonlinearImplicitSystem::ComputeResidual *	residual_object
void(*	jacobian )(const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)
NonlinearImplicitSystem::ComputeJacobian *	jacobian_object
void(*	matvec )(const NumericVector< Number > &X, NumericVector< Number > R, SparseMatrix< Number > J, sys_type &S)
NonlinearImplicitSystem::ComputeResidualandJacobian *	residual_and_jacobian_object
void(*	bounds )(NumericVector< Number > &XL, NumericVector< Number > &XU, sys_type &S)
NonlinearImplicitSystem::ComputeBounds *	bounds_object
void(*	nullspace )(std::vector< NumericVector< Number > * > &sp, sys_type &S)
NonlinearImplicitSystem::ComputeVectorSubspace *	nullspace_object
void(*	nearnullspace )(std::vector< NumericVector< Number > * > &sp, sys_type &S)
NonlinearImplicitSystem::ComputeVectorSubspace *	nearnullspace_object
void(*	user_presolve )(sys_type &S)
unsigned int	max_nonlinear_iterations
unsigned int	max_function_evaluations
Real	absolute_residual_tolerance
Real	relative_residual_tolerance
Real	absolute_step_tolerance
Real	relative_step_tolerance
unsigned int	max_linear_iterations
Real	initial_linear_tolerance
Real	minimum_linear_tolerance
bool	converged
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
void	build_mat_null_space (NonlinearImplicitSystem::ComputeVectorSubspace computeSubspaceObject, void()(std::vector< NumericVector< Number > * > &, sys_type &), MatNullSpace *)
void	increment_constructor_count (const std::string &name)
void	increment_destructor_count (const std::string &name)
Protected Attributes
SNES	_snes
SNESConvergedReason	_reason
PetscInt	_n_linear_iterations
unsigned	_current_nonlinear_iteration_number
bool	_zero_out_residual
bool	_zero_out_jacobian
bool	_default_monitor
sys_type &	_system
bool	_is_initialized
Preconditioner< T > *	_preconditioner
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
Friends
PetscErrorCode	__libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void *ctx)
PetscErrorCode	__libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat jac, Mat pc, MatStructure msflag, void ctx)
PetscErrorCode	__libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat jac, Mat pc, void *ctx)

Detailed Description

template<typename T>
class libMesh::PetscNonlinearSolver< T >

This class provides an interface to PETSc iterative solvers that is compatible with the libMesh NonlinearSolver<>

Author:: Benjamin Kirk, 2002-2007

Definition at line 63 of file petsc_nonlinear_solver.h.

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts [protected, inherited]

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

Definition at line 113 of file reference_counter.h.

template<typename T>

typedef NonlinearImplicitSystem libMesh::PetscNonlinearSolver< T >::sys_type

The type of system

Reimplemented from libMesh::NonlinearSolver< T >.

Definition at line 69 of file petsc_nonlinear_solver.h.

Constructor & Destructor Documentation

template<typename T>

libMesh::PetscNonlinearSolver< T >::PetscNonlinearSolver ( sys_type & system ) [explicit]

Constructor. Initializes Petsc data structures

template<typename T >

libMesh::PetscNonlinearSolver< T >::~PetscNonlinearSolver ( )

Destructor.

Definition at line 276 of file petsc_nonlinear_solver.C.

{
  this->clear ();
}

Member Function Documentation

template<typename T>

void libMesh::NonlinearSolver< T >::attach_preconditioner ( Preconditioner< T > * preconditioner ) [inherited]

Attaches a Preconditioner object to be used during the linear solves.

Definition at line 74 of file nonlinear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized.

{
  if (this->_is_initialized)
    libmesh_error_msg("Preconditioner must be attached before the solver is initialized!");

  _preconditioner = preconditioner;
}

template<typename T >

UniquePtr< NonlinearSolver< T > > libMesh::NonlinearSolver< T >::build	(	sys_type &	s,
		const SolverPackage	solver_package = `libMesh::default_solver_package()`
	)		`[static, inherited]`

Builds a NonlinearSolver using the nonlinear solver package specified by solver_package

Definition at line 37 of file nonlinear_solver.C.

References libMesh::PETSC_SOLVERS, and libMesh::TRILINOS_SOLVERS.

{
  // Build the appropriate solver
  switch (solver_package)
    {

#ifdef LIBMESH_HAVE_PETSC
    case PETSC_SOLVERS:
      return UniquePtr<NonlinearSolver<T> >(new PetscNonlinearSolver<T>(s));
#endif // LIBMESH_HAVE_PETSC

#ifdef LIBMESH_HAVE_NOX
    case TRILINOS_SOLVERS:
      return UniquePtr<NonlinearSolver<T> >(new NoxNonlinearSolver<T>(s));
#endif

    default:
      libmesh_error_msg("ERROR:  Unrecognized solver package: " << solver_package);
    }

  libmesh_error_msg("We'll never get here!");
  return UniquePtr<NonlinearSolver<T> >();
}

template<typename T >

void libMesh::PetscNonlinearSolver< T >::build_mat_null_space	(	NonlinearImplicitSystem::ComputeVectorSubspace *	computeSubspaceObject,
		void()(std::vector< NumericVector< Number > > &, sys_type &)	computeSubspace,
		MatNullSpace *	msp
	)		`[protected]`

Definition at line 396 of file petsc_nonlinear_solver.C.

References libMesh::ierr, and libMesh::PetscVector< T >::vec().

{
  PetscErrorCode ierr;
  std::vector<NumericVector<Number>* > sp;
  if (computeSubspaceObject)
    (*computeSubspaceObject)(sp, this->system());
  else
    (*computeSubspace)(sp, this->system());

  *msp = PETSC_NULL;
  if (sp.size())
    {
      Vec *modes;
      PetscScalar *dots;
      PetscInt nmodes = cast_int<PetscInt>(sp.size());

#if PETSC_RELEASE_LESS_THAN(3,5,0)
      ierr = PetscMalloc2(nmodes,Vec,&modes,nmodes,PetscScalar,&dots);
#else
      ierr = PetscMalloc2(nmodes,&modes,nmodes,&dots);
#endif
      LIBMESH_CHKERRABORT(ierr);

      for (PetscInt i=0; i<nmodes; ++i)
        {
          PetscVector<T>* pv = cast_ptr<PetscVector<T>*>(sp[i]);
          Vec v = pv->vec();

          ierr = VecDuplicate(v, modes+i);
          LIBMESH_CHKERRABORT(ierr);

          ierr = VecCopy(v,modes[i]);
          LIBMESH_CHKERRABORT(ierr);
        }

      // Normalize.
      ierr = VecNormalize(modes[0],PETSC_NULL);
      LIBMESH_CHKERRABORT(ierr);

      for (PetscInt i=1; i<nmodes; i++)
        {
          // Orthonormalize vec[i] against vec[0:i-1]
          ierr = VecMDot(modes[i],i,modes,dots);
          LIBMESH_CHKERRABORT(ierr);

          for (PetscInt j=0; j<i; j++)
            dots[j] *= -1.;

          ierr = VecMAXPY(modes[i],i,dots,modes);
          LIBMESH_CHKERRABORT(ierr);

          ierr = VecNormalize(modes[i],PETSC_NULL);
          LIBMESH_CHKERRABORT(ierr);
        }

      ierr = MatNullSpaceCreate(this->comm().get(), PETSC_FALSE, nmodes, modes, msp);
      LIBMESH_CHKERRABORT(ierr);

      for (PetscInt i=0; i<nmodes; ++i)
        {
          ierr = VecDestroy(modes+i);
          LIBMESH_CHKERRABORT(ierr);
        }

      ierr = PetscFree2(modes,dots);
      LIBMESH_CHKERRABORT(ierr);
    }
}

template<typename T >

void libMesh::PetscNonlinearSolver< T >::clear ( ) [virtual]

Release all memory and clear data structures.

Reimplemented from libMesh::NonlinearSolver< T >.

Definition at line 284 of file petsc_nonlinear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, and libMesh::initialized().

{
  if (this->initialized())
    {
      this->_is_initialized = false;

      PetscErrorCode ierr=0;

      ierr = LibMeshSNESDestroy(&_snes);
      LIBMESH_CHKERRABORT(ierr);

      // Reset the nonlinear iteration counter.  This information is only relevant
      // *during* the solve().  After the solve is completed it should return to
      // the default value of 0.
      _current_nonlinear_iteration_number = 0;
    }
}

const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const [inline, inherited]

Returns:: a reference to the Parallel::Communicator object used by this mesh.

Definition at line 86 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

  { return _communicator; }

void libMesh::ReferenceCounter::disable_print_counter_info ( ) [static, inherited]

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

{
  _enable_print_counter = false;
  return;
}

void libMesh::ReferenceCounter::enable_print_counter_info ( ) [static, inherited]

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

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

{
  _enable_print_counter = true;
  return;
}

template<typename T >

SNESConvergedReason libMesh::PetscNonlinearSolver< T >::get_converged_reason ( )

Returns the currently-available (or most recently obtained, if the SNES object has been destroyed) convergence reason. Refer to PETSc docs for the meaning of different SNESConvergedReasons.

Definition at line 627 of file petsc_nonlinear_solver.C.

References libMesh::initialized().

{
  PetscErrorCode ierr=0;

  if (this->initialized())
    {
      ierr = SNESGetConvergedReason(_snes, &_reason);
      LIBMESH_CHKERRABORT(ierr);
    }

  return _reason;
}

template<typename T>

virtual unsigned libMesh::PetscNonlinearSolver< T >::get_current_nonlinear_iteration_number ( ) const [inline, virtual]

If called *during* the solve(), for example by the user-specified residual or Jacobian function, returns the current nonlinear iteration number.

Implements libMesh::NonlinearSolver< T >.

Definition at line 131 of file petsc_nonlinear_solver.h.

References libMesh::PetscNonlinearSolver< T >::_current_nonlinear_iteration_number.

{ return _current_nonlinear_iteration_number; }

std::string libMesh::ReferenceCounter::get_info ( ) [static, inherited]

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().

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

  std::ostringstream oss;

  oss << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";

  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;

      oss << "| " << name << " reference count information:\n"
          << "|  Creations:    " << creations    << '\n'
          << "|  Destructions: " << destructions << '\n';
    }

  oss << " ---------------------------------------------------------------------------- \n";

  return oss.str();

#else

  return "";

#endif
}

template<typename T >

int libMesh::PetscNonlinearSolver< T >::get_total_linear_iterations ( ) [virtual]

Get the total number of linear iterations done in the last solve

Implements libMesh::NonlinearSolver< T >.

Definition at line 641 of file petsc_nonlinear_solver.C.

{
  return _n_linear_iterations;
}

void libMesh::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 163 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

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

  p.first++;
}

void libMesh::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 176 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

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

  p.second++;
}

template<typename T >

void libMesh::PetscNonlinearSolver< T >::init ( const char * name = NULL ) [virtual]

Initialize data structures if not done so already. May assign a name to the solver in some implementations

Implements libMesh::NonlinearSolver< T >.

Definition at line 305 of file petsc_nonlinear_solver.C.

References libMesh::__libmesh_petsc_preconditioner_apply(), libMesh::__libmesh_petsc_preconditioner_setup(), libMesh::__libmesh_petsc_snes_monitor(), libMesh::libMeshPrivateData::_is_initialized, DMlibMeshSetSystem(), libMesh::initialized(), and libMesh::pc.

Referenced by libMesh::PetscNonlinearSolver< T >::snes().

{
  // Initialize the data structures if not done so already.
  if (!this->initialized())
    {
      this->_is_initialized = true;

      PetscErrorCode ierr=0;

#if PETSC_VERSION_LESS_THAN(2,1,2)
      // At least until Petsc 2.1.1, the SNESCreate had a different calling syntax.
      // The second argument was of type SNESProblemType, and could have a value of
      // either SNES_NONLINEAR_EQUATIONS or SNES_UNCONSTRAINED_MINIMIZATION.
      ierr = SNESCreate(this->comm().get(), SNES_NONLINEAR_EQUATIONS, &_snes);
      LIBMESH_CHKERRABORT(ierr);

#else

      ierr = SNESCreate(this->comm().get(),&_snes);
      LIBMESH_CHKERRABORT(ierr);

#endif

      if (name)
        {
          ierr = SNESSetOptionsPrefix(_snes, name);
          LIBMESH_CHKERRABORT(ierr);
        }

#if !PETSC_RELEASE_LESS_THAN(3,3,0)
      // Attaching a DM to SNES.
      DM dm;
      ierr = DMCreate(this->comm().get(), &dm);LIBMESH_CHKERRABORT(ierr);
      ierr = DMSetType(dm,DMLIBMESH);LIBMESH_CHKERRABORT(ierr);
      ierr = DMlibMeshSetSystem(dm,this->system());LIBMESH_CHKERRABORT(ierr);
      if (name)
        {
          ierr = DMSetOptionsPrefix(dm,name);    LIBMESH_CHKERRABORT(ierr);
        }
      ierr = DMSetFromOptions(dm);               LIBMESH_CHKERRABORT(ierr);
      ierr = DMSetUp(dm);                        LIBMESH_CHKERRABORT(ierr);
      ierr = SNESSetDM(this->_snes, dm);         LIBMESH_CHKERRABORT(ierr);
      // SNES now owns the reference to dm.
      ierr = DMDestroy(&dm);                     LIBMESH_CHKERRABORT(ierr);

#endif

      if (_default_monitor)
        {
#if PETSC_VERSION_LESS_THAN(2,3,3)
          ierr = SNESSetMonitor (_snes, __libmesh_petsc_snes_monitor,
                                 this, PETSC_NULL);
#else
          // API name change in PETSc 2.3.3
          ierr = SNESMonitorSet (_snes, __libmesh_petsc_snes_monitor,
                                 this, PETSC_NULL);
#endif
          LIBMESH_CHKERRABORT(ierr);
        }

#if PETSC_VERSION_LESS_THAN(3,1,0)
      // Cannot call SNESSetOptions before SNESSetFunction when using
      // any matrix free options with PETSc 3.1.0+
      ierr = SNESSetFromOptions(_snes);
      LIBMESH_CHKERRABORT(ierr);
#endif

      if(this->_preconditioner)
        {
          KSP ksp;
          ierr = SNESGetKSP (_snes, &ksp);
          LIBMESH_CHKERRABORT(ierr);
          PC pc;
          ierr = KSPGetPC(ksp,&pc);
          LIBMESH_CHKERRABORT(ierr);

          this->_preconditioner->init();

          PCSetType(pc, PCSHELL);
          PCShellSetContext(pc,(void*)this->_preconditioner);

          //Re-Use the shell functions from petsc_linear_solver
          PCShellSetSetUp(pc,__libmesh_petsc_preconditioner_setup);
          PCShellSetApply(pc,__libmesh_petsc_preconditioner_apply);
        }
    }
}

template<typename T>

bool libMesh::NonlinearSolver< T >::initialized ( ) const [inline, inherited]

Returns:: true if the data structures are initialized, false otherwise.

Definition at line 82 of file nonlinear_solver.h.

{ return _is_initialized; }

static unsigned int libMesh::ReferenceCounter::n_objects ( ) [inline, static, inherited]

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

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

  { return _n_objects; }

processor_id_type libMesh::ParallelObject::n_processors ( ) const [inline, inherited]

Returns:: the number of processors in the group.

Definition at line 92 of file parallel_object.h.

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

  { return cast_int<processor_id_type>(_communicator.size()); }

template<typename T >

void libMesh::PetscNonlinearSolver< T >::print_converged_reason ( ) [virtual]

Prints a useful message about why the latest nonlinear solve con(di)verged.

Reimplemented from libMesh::NonlinearSolver< T >.

Definition at line 617 of file petsc_nonlinear_solver.C.

References libMesh::out.

{

  libMesh::out << "Nonlinear solver convergence/divergence reason: "
               << SNESConvergedReasons[this->get_converged_reason()] << std::endl;
}

void libMesh::ReferenceCounter::print_info ( std::ostream & out = libMesh::out ) [static, inherited]

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().

{
  if( _enable_print_counter ) out_stream << ReferenceCounter::get_info();
}

processor_id_type libMesh::ParallelObject::processor_id ( ) const [inline, inherited]

Returns:: the rank of this processor in the group.

Definition at line 98 of file parallel_object.h.

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

  { return cast_int<processor_id_type>(_communicator.rank()); }

template<typename T>

void libMesh::PetscNonlinearSolver< T >::set_jacobian_zero_out ( bool state ) [inline]

Set if the jacobian should be zeroed out in the callback

Definition at line 141 of file petsc_nonlinear_solver.h.

References libMesh::PetscNonlinearSolver< T >::_zero_out_jacobian.

{ _zero_out_jacobian = state; }

template<typename T>

void libMesh::PetscNonlinearSolver< T >::set_residual_zero_out ( bool state ) [inline]

Set if the residual should be zeroed out in the callback

Definition at line 136 of file petsc_nonlinear_solver.h.

References libMesh::PetscNonlinearSolver< T >::_zero_out_residual.

{ _zero_out_residual = state; }

template<typename T>

SNES libMesh::PetscNonlinearSolver< T >::snes ( ) [inline]

Returns the raw PETSc snes context pointer.

Definition at line 96 of file petsc_nonlinear_solver.h.

References libMesh::PetscNonlinearSolver< T >::_snes, and libMesh::PetscNonlinearSolver< T >::init().

{ this->init(); return _snes; }

template<typename T >

std::pair< unsigned int, Real > libMesh::PetscNonlinearSolver< T >::solve	(	SparseMatrix< T > &	jac_in,
		NumericVector< T > &	x_in,
		NumericVector< T > &	r_in,
		const double	,
		const unsigned	int
	)		`[virtual]`

Call the Petsc solver. It calls the method below, using the same matrix for the system and preconditioner matrices.

Implements libMesh::NonlinearSolver< T >.

Definition at line 470 of file petsc_nonlinear_solver.C.

References libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_residual(), libMesh::TriangleWrapper::init(), libMesh::PetscMatrix< T >::mat(), libMesh::Real, and libMesh::START_LOG().

{
  START_LOG("solve()", "PetscNonlinearSolver");
  this->init ();

  // Make sure the data passed in are really of Petsc types
  PetscMatrix<T>* jac = cast_ptr<PetscMatrix<T>*>(&jac_in);
  PetscVector<T>* x   = cast_ptr<PetscVector<T>*>(&x_in);
  PetscVector<T>* r   = cast_ptr<PetscVector<T>*>(&r_in);

  PetscErrorCode ierr=0;
  PetscInt n_iterations =0;
  // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal
  Real final_residual_norm=0.;

  ierr = SNESSetFunction (_snes, r->vec(), __libmesh_petsc_snes_residual, this);
  LIBMESH_CHKERRABORT(ierr);

  // Only set the jacobian function if we've been provided with something to call.
  // This allows a user to set their own jacobian function if they want to
  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
    {
      ierr = SNESSetJacobian (_snes, jac->mat(), jac->mat(), __libmesh_petsc_snes_jacobian, this);
      LIBMESH_CHKERRABORT(ierr);
    }
#if !PETSC_VERSION_LESS_THAN(3,3,0)
  // Only set the nullspace if we have a way of computing it and the result is non-empty.
  if (this->nullspace || this->nullspace_object)
    {
      MatNullSpace msp;
      this->build_mat_null_space(this->nullspace_object, this->nullspace, &msp);
      if (msp)
        {
          ierr = MatSetNullSpace(jac->mat(), msp);
          LIBMESH_CHKERRABORT(ierr);

          ierr = MatNullSpaceDestroy(&msp);
          LIBMESH_CHKERRABORT(ierr);
        }
    }

  // Only set the nearnullspace if we have a way of computing it and the result is non-empty.
  if (this->nearnullspace || this->nearnullspace_object)
    {
      MatNullSpace msp = PETSC_NULL;
      this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, &msp);

      if(msp) {
        ierr = MatSetNearNullSpace(jac->mat(), msp);
        LIBMESH_CHKERRABORT(ierr);

        ierr = MatNullSpaceDestroy(&msp);
        LIBMESH_CHKERRABORT(ierr);
      }
    }
#endif
  // Have the Krylov subspace method use our good initial guess rather than 0
  KSP ksp;
  ierr = SNESGetKSP (_snes, &ksp);
  LIBMESH_CHKERRABORT(ierr);

  // Set the tolerances for the iterative solver.  Use the user-supplied
  // tolerance for the relative residual & leave the others at default values
  ierr = KSPSetTolerances (ksp, this->initial_linear_tolerance, PETSC_DEFAULT,
                           PETSC_DEFAULT, this->max_linear_iterations);
  LIBMESH_CHKERRABORT(ierr);

  // Set the tolerances for the non-linear solver.
  ierr = SNESSetTolerances(_snes, this->absolute_residual_tolerance, this->relative_residual_tolerance,
                           this->relative_step_tolerance, this->max_nonlinear_iterations, this->max_function_evaluations);
  LIBMESH_CHKERRABORT(ierr);

  //Pull in command-line options
  KSPSetFromOptions(ksp);
  SNESSetFromOptions(_snes);

  if (this->user_presolve)
    this->user_presolve(this->system());

  //Set the preconditioning matrix
  if(this->_preconditioner)
    {
      this->_preconditioner->set_matrix(jac_in);
      this->_preconditioner->init();
    }

  //    ierr = KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
  //           LIBMESH_CHKERRABORT(ierr);

  // Older versions (at least up to 2.1.5) of SNESSolve took 3 arguments,
  // the last one being a pointer to an int to hold the number of iterations required.
# if PETSC_VERSION_LESS_THAN(2,2,0)

  ierr = SNESSolve (_snes, x->vec(), &n_iterations);
  LIBMESH_CHKERRABORT(ierr);

  // 2.2.x style
#elif PETSC_VERSION_LESS_THAN(2,3,0)

  ierr = SNESSolve (_snes, x->vec());
  LIBMESH_CHKERRABORT(ierr);

  // 2.3.x & newer style
#else

  ierr = SNESSolve (_snes, PETSC_NULL, x->vec());
  LIBMESH_CHKERRABORT(ierr);

  ierr = SNESGetIterationNumber(_snes,&n_iterations);
  LIBMESH_CHKERRABORT(ierr);

  ierr = SNESGetLinearSolveIterations(_snes, &_n_linear_iterations);
  LIBMESH_CHKERRABORT(ierr);

#endif

  // SNESGetFunction has been around forever and should work on all
  // versions of PETSc.  This is also now the recommended approach
  // according to the documentation for the PETSc 3.5.1 release:
  // http://www.mcs.anl.gov/petsc/documentation/changes/35.html
  Vec f;
  ierr = SNESGetFunction(_snes, &f, 0, 0);
  LIBMESH_CHKERRABORT(ierr);
  ierr = VecNorm(f, NORM_2, &final_residual_norm);
  LIBMESH_CHKERRABORT(ierr);

  // Get and store the reason for convergence
  SNESGetConvergedReason(_snes, &_reason);

  //Based on Petsc 2.3.3 documentation all diverged reasons are negative
  this->converged = (_reason >= 0);

  this->clear();

  STOP_LOG("solve()", "PetscNonlinearSolver");

  // return the # of its. and the final residual norm.
  return std::make_pair(n_iterations, final_residual_norm);
}

template<typename T>

const sys_type& libMesh::NonlinearSolver< T >::system ( ) const [inline, inherited]

Returns:: a constant reference to the system we are solving.

Definition at line 216 of file nonlinear_solver.h.

Referenced by libMesh::__libmesh_petsc_snes_residual(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), and libMesh::Problem_Interface::computePreconditioner().

{ return _system; }

template<typename T>

sys_type& libMesh::NonlinearSolver< T >::system ( ) [inline, inherited]

Returns:: a writeable reference to the system we are solving.

Definition at line 221 of file nonlinear_solver.h.

{ return _system; }

template<typename T>

void libMesh::PetscNonlinearSolver< T >::use_default_monitor ( bool state ) [inline]

Set to true to use the libMeash's default monitor, set to false to use your own

Definition at line 146 of file petsc_nonlinear_solver.h.

References libMesh::PetscNonlinearSolver< T >::_default_monitor.

{ _default_monitor = state; }

Friends And Related Function Documentation

template<typename T>

PetscErrorCode __libmesh_petsc_snes_jacobian	(	SNES	snes,
		Vec	x,
		Mat *	jac,
		Mat *	pc,
		MatStructure *	msflag,
		void *	ctx
	)		`[friend]`

template<typename T>

PetscErrorCode __libmesh_petsc_snes_jacobian	(	SNES	snes,
		Vec	x,
		Mat	jac,
		Mat	pc,
		void *	ctx
	)		`[friend]`

template<typename T>

PetscErrorCode __libmesh_petsc_snes_residual	(	SNES	snes,
		Vec	x,
		Vec	r,
		void *	ctx
	)		`[friend]`

Definition at line 81 of file petsc_nonlinear_solver.C.

  {
    START_LOG("residual()", "PetscNonlinearSolver");

    PetscErrorCode ierr=0;

    libmesh_assert(x);
    libmesh_assert(r);
    libmesh_assert(ctx);

    // No way to safety-check this cast, since we got a void*...
    PetscNonlinearSolver<Number>* solver =
      static_cast<PetscNonlinearSolver<Number>*> (ctx);

    // Get the current iteration number from the snes object,
    // store it in the PetscNonlinearSolver object for possible use
    // by the user's residual function.
    {
      PetscInt n_iterations = 0;
      ierr = SNESGetIterationNumber(snes, &n_iterations);
      CHKERRABORT(solver->comm().get(),ierr);
      solver->_current_nonlinear_iteration_number = cast_int<unsigned>(n_iterations);
    }

    NonlinearImplicitSystem &sys = solver->system();

    PetscVector<Number>& X_sys = *cast_ptr<PetscVector<Number>*>(sys.solution.get());
    PetscVector<Number>& R_sys = *cast_ptr<PetscVector<Number>*>(sys.rhs);
    PetscVector<Number> X_global(x, sys.comm()), R(r, sys.comm());

    // Use the systems update() to get a good local version of the parallel solution
    X_global.swap(X_sys);
    R.swap(R_sys);

    sys.get_dof_map().enforce_constraints_exactly(sys);

    sys.update();

    //Swap back
    X_global.swap(X_sys);
    R.swap(R_sys);

    if (solver->_zero_out_residual)
      R.zero();

    //-----------------------------------------------------------------------------
    // if the user has provided both function pointers and objects only the pointer
    // will be used, so catch that as an error
    if (solver->residual && solver->residual_object)
      libmesh_error_msg("ERROR: cannot specifiy both a function and object to compute the Residual!");

    if (solver->matvec && solver->residual_and_jacobian_object)
      libmesh_error_msg("ERROR: cannot specifiy both a function and object to compute the combined Residual & Jacobian!");

    if (solver->residual != NULL)
      solver->residual(*sys.current_local_solution.get(), R, sys);

    else if (solver->residual_object != NULL)
      solver->residual_object->residual(*sys.current_local_solution.get(), R, sys);

    else if (solver->matvec != NULL)
      solver->matvec (*sys.current_local_solution.get(), &R, NULL, sys);

    else if (solver->residual_and_jacobian_object != NULL)
      solver->residual_and_jacobian_object->residual_and_jacobian (*sys.current_local_solution.get(), &R, NULL, sys);

    else
      libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");

    R.close();
    X_global.close();

    STOP_LOG("residual()", "PetscNonlinearSolver");

    return ierr;
  }

Member Data Documentation

const Parallel::Communicator& libMesh::ParallelObject::_communicator [protected, inherited]

Definition at line 104 of file parallel_object.h.

Referenced by libMesh::EquationSystems::build_solution_vector(), libMesh::ParallelObject::comm(), libMesh::EquationSystems::get_solution(), libMesh::ParallelObject::n_processors(), libMesh::ParallelObject::operator=(), and libMesh::ParallelObject::processor_id().

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts [static, protected, inherited]

Actually holds the data.

Definition at line 118 of file reference_counter.h.

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

template<typename T>

unsigned libMesh::PetscNonlinearSolver< T >::_current_nonlinear_iteration_number [protected]

Stores the current nonlinear iteration number

Definition at line 171 of file petsc_nonlinear_solver.h.

Referenced by libMesh::__libmesh_petsc_snes_residual(), and libMesh::PetscNonlinearSolver< T >::get_current_nonlinear_iteration_number().

template<typename T>

bool libMesh::PetscNonlinearSolver< T >::_default_monitor [protected]

true if we want the default monitor to be set, false for no monitor (i.e. user code can use their own)

Definition at line 186 of file petsc_nonlinear_solver.h.

Referenced by libMesh::PetscNonlinearSolver< T >::use_default_monitor().

bool libMesh::ReferenceCounter::_enable_print_counter = true [static, protected, inherited]

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

Definition at line 137 of file reference_counter.h.

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

template<typename T>

bool libMesh::NonlinearSolver< T >::_is_initialized [protected, inherited]

Flag indicating if the data structures have been initialized.

Definition at line 297 of file nonlinear_solver.h.

Referenced by libMesh::NonlinearSolver< Number >::initialized().

Threads::spin_mutex libMesh::ReferenceCounter::_mutex [static, protected, inherited]

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 131 of file reference_counter.h.

template<typename T>

PetscInt libMesh::PetscNonlinearSolver< T >::_n_linear_iterations [protected]

Stores the total number of linear iterations from the last solve.

Definition at line 166 of file petsc_nonlinear_solver.h.

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects [static, protected, inherited]

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

Definition at line 126 of file reference_counter.h.

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

template<typename T>

Preconditioner<T>* libMesh::NonlinearSolver< T >::_preconditioner [protected, inherited]

Holds the Preconditioner object to be used for the linear solves.

Definition at line 302 of file nonlinear_solver.h.

template<typename T>

SNESConvergedReason libMesh::PetscNonlinearSolver< T >::_reason [protected]

Store the reason for SNES convergence/divergence for use even after the _snes has been cleared. Note that print_converged_reason() will always *try* to get the current reason with SNESGetConvergedReason(), but if the SNES object has already been cleared, it will fall back on this stored value. Note that this value is therefore necessarily *not* cleared by the clear() function.

Definition at line 161 of file petsc_nonlinear_solver.h.

template<typename T>

SNES libMesh::PetscNonlinearSolver< T >::_snes [protected]

Nonlinear solver context

Definition at line 152 of file petsc_nonlinear_solver.h.

Referenced by libMesh::PetscNonlinearSolver< T >::snes().

template<typename T>

sys_type& libMesh::NonlinearSolver< T >::_system [protected, inherited]

A reference to the system we are solving.

Definition at line 292 of file nonlinear_solver.h.

Referenced by libMesh::NonlinearSolver< Number >::system().

template<typename T>

bool libMesh::PetscNonlinearSolver< T >::_zero_out_jacobian [protected]

true to zero out jacobian before going into application level call-back, otherwise false

Definition at line 181 of file petsc_nonlinear_solver.h.

Referenced by libMesh::PetscNonlinearSolver< T >::set_jacobian_zero_out().

template<typename T>

bool libMesh::PetscNonlinearSolver< T >::_zero_out_residual [protected]

true to zero out residual before going into application level call-back, otherwise false

Definition at line 176 of file petsc_nonlinear_solver.h.

Referenced by libMesh::__libmesh_petsc_snes_residual(), and libMesh::PetscNonlinearSolver< T >::set_residual_zero_out().

template<typename T>

Real libMesh::NonlinearSolver< T >::absolute_residual_tolerance [inherited]

The NonlinearSolver 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 248 of file nonlinear_solver.h.

template<typename T>

Real libMesh::NonlinearSolver< T >::absolute_step_tolerance [inherited]

The NonlinearSolver 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.

Note that not all NonlinearSolvers support relative_step_tolerance!

Definition at line 262 of file nonlinear_solver.h.

template<typename T>

void(* libMesh::NonlinearSolver< T >::bounds)(NumericVector< Number > &XL, NumericVector< Number > &XU, sys_type &S) [inherited]

Function that computes the lower and upper bounds XL and XU on the solution of the nonlinear system.

Definition at line 172 of file nonlinear_solver.h.

template<typename T>

NonlinearImplicitSystem::ComputeBounds* libMesh::NonlinearSolver< T >::bounds_object [inherited]

Object that computes the bounds vectors $ XL $ and $ XU $ .

Definition at line 178 of file nonlinear_solver.h.

template<typename T>

bool libMesh::NonlinearSolver< T >::converged [inherited]

After a call to solve this will reflect whether or not the nonlinear solve was successful.

Definition at line 286 of file nonlinear_solver.h.

template<typename T>

Real libMesh::NonlinearSolver< T >::initial_linear_tolerance [inherited]

Any required linear solves will at first be done with this tolerance; the NonlinearSolver may tighten the tolerance for later solves.

Definition at line 275 of file nonlinear_solver.h.

template<typename T>

void(* libMesh::NonlinearSolver< T >::jacobian)(const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S) [inherited]

Function that computes the Jacobian J(X) of the nonlinear system at the input iterate X.

Definition at line 140 of file nonlinear_solver.h.

Referenced by libMesh::Problem_Interface::computeJacobian(), and libMesh::Problem_Interface::computePreconditioner().

template<typename T>

NonlinearImplicitSystem::ComputeJacobian* libMesh::NonlinearSolver< T >::jacobian_object [inherited]

Object that computes the Jacobian J(X) of the nonlinear system at the input iterate X.

Definition at line 148 of file nonlinear_solver.h.

Referenced by libMesh::Problem_Interface::computeJacobian(), and libMesh::Problem_Interface::computePreconditioner().

template<typename T>

void(* libMesh::NonlinearSolver< T >::matvec)(const NumericVector< Number > &X, NumericVector< Number > *R, SparseMatrix< Number > *J, sys_type &S) [inherited]

Function that computes either the residual $ R(X) $ or the Jacobian $ J(X) $ of the nonlinear system at the input iterate $ X $ . Note that either R or J could be XSNULL.

Definition at line 156 of file nonlinear_solver.h.

Referenced by libMesh::__libmesh_petsc_snes_residual(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), and libMesh::Problem_Interface::computePreconditioner().

template<typename T>

unsigned int libMesh::NonlinearSolver< T >::max_function_evaluations [inherited]

Maximum number of function evaluations.

Definition at line 236 of file nonlinear_solver.h.

template<typename T>

unsigned int libMesh::NonlinearSolver< T >::max_linear_iterations [inherited]

Each linear solver step should exit after max_linear_iterations is exceeded.

Definition at line 269 of file nonlinear_solver.h.

template<typename T>

unsigned int libMesh::NonlinearSolver< T >::max_nonlinear_iterations [inherited]

Maximum number of non-linear iterations.

Definition at line 231 of file nonlinear_solver.h.

template<typename T>

Real libMesh::NonlinearSolver< T >::minimum_linear_tolerance [inherited]

The tolerance for linear solves is kept above this minimum

Definition at line 280 of file nonlinear_solver.h.

template<typename T>

void(* libMesh::NonlinearSolver< T >::nearnullspace)(std::vector< NumericVector< Number > * > &sp, sys_type &S) [inherited]

Function that computes a basis for the Jacobian's near nullspace -- the set of "low energy modes" -- that can be used for AMG coarsening, if the solver supports it (e.g., ML, PETSc's GAMG).

Definition at line 201 of file nonlinear_solver.h.

template<typename T>

NonlinearImplicitSystem::ComputeVectorSubspace* libMesh::NonlinearSolver< T >::nearnullspace_object [inherited]

A callable object that computes a basis for the Jacobian's near nullspace -- the set of "low energy modes" -- that can be used for AMG coarsening, if the solver supports it (e.g., ML, PETSc's GAMG).

Definition at line 208 of file nonlinear_solver.h.

template<typename T>

void(* libMesh::NonlinearSolver< T >::nullspace)(std::vector< NumericVector< Number > * > &sp, sys_type &S) [inherited]

Function that computes a basis for the Jacobian's nullspace -- the kernel or the "zero energy modes" -- that can be used in solving a degenerate problem iteratively, if the solver supports it (e.g., PETSc's KSP).

Definition at line 186 of file nonlinear_solver.h.

template<typename T>

NonlinearImplicitSystem::ComputeVectorSubspace* libMesh::NonlinearSolver< T >::nullspace_object [inherited]

A callable object that computes a basis for the Jacobian's nullspace -- the kernel or the "zero energy modes" -- that can be used in solving a degenerate problem iteratively, if the solver supports it (e.g., PETSc's KSP).

Definition at line 194 of file nonlinear_solver.h.

template<typename T>

Real libMesh::NonlinearSolver< T >::relative_residual_tolerance [inherited]

Definition at line 249 of file nonlinear_solver.h.

template<typename T>

Real libMesh::NonlinearSolver< T >::relative_step_tolerance [inherited]

Definition at line 263 of file nonlinear_solver.h.

template<typename T>

void(* libMesh::NonlinearSolver< T >::residual)(const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S) [inherited]

Function that computes the residual R(X) of the nonlinear system at the input iterate X.

Definition at line 126 of file nonlinear_solver.h.

Referenced by libMesh::__libmesh_petsc_snes_residual(), and libMesh::Problem_Interface::computeF().

template<typename T>

NonlinearImplicitSystem::ComputeResidualandJacobian* libMesh::NonlinearSolver< T >::residual_and_jacobian_object [inherited]

Object that computes either the residual $ R(X) $ or the Jacobian $ J(X) $ of the nonlinear system at the input iterate $ X $ . Note that either R or J could be XSNULL.

Definition at line 167 of file nonlinear_solver.h.

Referenced by libMesh::__libmesh_petsc_snes_residual(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), and libMesh::Problem_Interface::computePreconditioner().

template<typename T>

NonlinearImplicitSystem::ComputeResidual* libMesh::NonlinearSolver< T >::residual_object [inherited]

Object that computes the residual R(X) of the nonlinear system at the input iterate X.

Definition at line 134 of file nonlinear_solver.h.

Referenced by libMesh::__libmesh_petsc_snes_residual(), and libMesh::Problem_Interface::computeF().

template<typename T>

void(* libMesh::NonlinearSolver< T >::user_presolve)(sys_type &S) [inherited]

Definition at line 211 of file nonlinear_solver.h.

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