trilinos_nox_nonlinear_solver.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#include "libmesh/libmesh_common.h"

#if defined(LIBMESH_TRILINOS_HAVE_NOX) && defined(LIBMESH_TRILINOS_HAVE_EPETRA)

// Local Includes
#include "libmesh/libmesh_logging.h"
#include "libmesh/dof_map.h"
#include "libmesh/nonlinear_implicit_system.h"
#include "libmesh/trilinos_nox_nonlinear_solver.h"
#include "libmesh/system.h"
#include "libmesh/trilinos_epetra_vector.h"
#include "libmesh/trilinos_epetra_matrix.h"
#include "libmesh/trilinos_preconditioner.h"

// Trilinos includes
#include "libmesh/ignore_warnings.h"
#include "NOX_Epetra_MatrixFree.H"
#include "NOX_Epetra_LinearSystem_AztecOO.H"
#include "NOX_Epetra_Group.H"
#include "NOX_Epetra_Vector.H"
#include "libmesh/restore_warnings.h"

namespace libMesh
{

class Problem_Interface : public NOX::Epetra::Interface::Required,
                          public NOX::Epetra::Interface::Jacobian,
                          public NOX::Epetra::Interface::Preconditioner
{
public:
  explicit
  Problem_Interface(NoxNonlinearSolver<Number> * solver);

  ~Problem_Interface();

  bool computeF(const Epetra_Vector & x,
                Epetra_Vector & FVec,
                NOX::Epetra::Interface::Required::FillType fillType);

  bool computeJacobian(const Epetra_Vector & x,
                       Epetra_Operator & Jac);

  bool computePrecMatrix(const Epetra_Vector & x,
                         Epetra_RowMatrix & M);

  bool computePreconditioner(const Epetra_Vector & x,
                             Epetra_Operator & Prec,
                             Teuchos::ParameterList * p);

  NoxNonlinearSolver<Number> * _solver;
};



Problem_Interface::Problem_Interface(NoxNonlinearSolver<Number> * solver) :
  _solver(solver)
{}



Problem_Interface::~Problem_Interface()
{}



bool Problem_Interface::computeF(const Epetra_Vector & x,
                                 Epetra_Vector & r,
                                 NOX::Epetra::Interface::Required::FillType /*fillType*/)
{
  LOG_SCOPE("residual()", "TrilinosNoxNonlinearSolver");

  NonlinearImplicitSystem & sys = _solver->system();

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

  // 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);

  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 specify both a function and object to compute the Residual!");

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

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

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

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

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

  else
    return false;

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

  return true;
}



bool Problem_Interface::computeJacobian(const Epetra_Vector & x,
                                        Epetra_Operator & jac)
{
  LOG_SCOPE("jacobian()", "TrilinosNoxNonlinearSolver");

  NonlinearImplicitSystem & sys = _solver->system();

  EpetraMatrix<Number> J(&dynamic_cast<Epetra_FECrsMatrix &>(jac), sys.comm());
  EpetraVector<Number> & X_sys = *cast_ptr<EpetraVector<Number> *>(sys.solution.get());
  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm());

  // Set the dof maps
  J.attach_dof_map(sys.get_dof_map());

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

  sys.get_dof_map().enforce_constraints_exactly(sys);
  sys.update();

  X_global.swap(X_sys);

  //-----------------------------------------------------------------------------
  // if the user has provided both function pointers and objects only the pointer
  // will be used, so catch that as an error
  if (_solver->jacobian && _solver->jacobian_object)
    libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Jacobian!");

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

  if (_solver->jacobian != nullptr)
    _solver->jacobian(*sys.current_local_solution.get(), J, sys);

  else if (_solver->jacobian_object != nullptr)
    _solver->jacobian_object->jacobian(*sys.current_local_solution.get(), J, sys);

  else if (_solver->matvec != nullptr)
    _solver->matvec(*sys.current_local_solution.get(), nullptr, &J, sys);

  else if (_solver->residual_and_jacobian_object != nullptr)
    _solver->residual_and_jacobian_object->residual_and_jacobian (*sys.current_local_solution.get(), nullptr, &J, sys);

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

  J.close();
  X_global.close();

  return true;
}



bool Problem_Interface::computePrecMatrix(const Epetra_Vector & /*x*/, Epetra_RowMatrix & /*M*/)
{
  //   libMesh::out << "ERROR: Problem_Interface::preconditionVector() - Use Explicit Jacobian only for this test problem!" << endl;
  throw 1;
}



bool Problem_Interface::computePreconditioner(const Epetra_Vector & x,
                                              Epetra_Operator & prec,
                                              Teuchos::ParameterList * /*p*/)
{
  LOG_SCOPE("preconditioner()", "TrilinosNoxNonlinearSolver");

  NonlinearImplicitSystem & sys = _solver->system();
  TrilinosPreconditioner<Number> & tpc = dynamic_cast<TrilinosPreconditioner<Number> &>(prec);

  EpetraMatrix<Number> J(dynamic_cast<Epetra_FECrsMatrix *>(tpc.mat()),sys.comm());
  EpetraVector<Number> & X_sys = *cast_ptr<EpetraVector<Number> *>(sys.solution.get());
  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm());

  // Set the dof maps
  J.attach_dof_map(sys.get_dof_map());

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

  sys.get_dof_map().enforce_constraints_exactly(sys);
  sys.update();

  X_global.swap(X_sys);

  //-----------------------------------------------------------------------------
  // if the user has provided both function pointers and objects only the pointer
  // will be used, so catch that as an error
  if (_solver->jacobian && _solver->jacobian_object)
    libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Jacobian!");

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

  if (_solver->jacobian != nullptr)
    _solver->jacobian(*sys.current_local_solution.get(), J, sys);

  else if (_solver->jacobian_object != nullptr)
    _solver->jacobian_object->jacobian(*sys.current_local_solution.get(), J, sys);

  else if (_solver->matvec != nullptr)
    _solver->matvec(*sys.current_local_solution.get(), nullptr, &J, sys);

  else if (_solver->residual_and_jacobian_object != nullptr)
    _solver->residual_and_jacobian_object->residual_and_jacobian (*sys.current_local_solution.get(), nullptr, &J, sys);

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

  J.close();
  X_global.close();

  tpc.compute();

  return true;
}



//---------------------------------------------------------------------
// NoxNonlinearSolver<> methods
template <typename T>
void NoxNonlinearSolver<T>::clear ()
{
}



template <typename T>
void NoxNonlinearSolver<T>::init (const char * /*name*/)
{
  if (!this->initialized())
    _interface = new Problem_Interface(this);
}



template <typename T>
std::pair<unsigned int, Real>
NoxNonlinearSolver<T>::solve (SparseMatrix<T> &  /* jac_in */,  // System Jacobian Matrix
                              NumericVector<T> & x_in,          // Solution vector
                              NumericVector<T> & /* r_in */,    // Residual vector
                              const double,                    // Stopping tolerance
                              const unsigned int)
{
  this->init ();

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

  EpetraVector<T> * x_epetra = cast_ptr<EpetraVector<T> *>(&x_in);
  // Creating a Teuchos::RCP as they do in NOX examples does not work here - we get some invalid memory references
  // thus we make a local copy
  NOX::Epetra::Vector x(*x_epetra->vec());

  Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr = Teuchos::rcp(new Teuchos::ParameterList);
  Teuchos::ParameterList & nlParams = *(nlParamsPtr.get());
  nlParams.set("Nonlinear Solver", "Line Search Based");

  //print params
  Teuchos::ParameterList & printParams = nlParams.sublist("Printing");
  printParams.set("Output Precision", 3);
  printParams.set("Output Processor", 0);
  printParams.set("Output Information",
                  NOX::Utils::OuterIteration +
                  NOX::Utils::OuterIterationStatusTest +
                  NOX::Utils::InnerIteration +
                  NOX::Utils::LinearSolverDetails +
                  NOX::Utils::Parameters +
                  NOX::Utils::Details +
                  NOX::Utils::Warning);

  Teuchos::ParameterList & dirParams = nlParams.sublist("Direction");
  dirParams.set("Method", "Newton");
  Teuchos::ParameterList & newtonParams = dirParams.sublist("Newton");
  newtonParams.set("Forcing Term Method", "Constant");

  Teuchos::ParameterList & lsParams = newtonParams.sublist("Linear Solver");
  lsParams.set("Aztec Solver", "GMRES");
  lsParams.set("Max Iterations", static_cast<int>(this->max_linear_iterations));
  lsParams.set("Tolerance", this->initial_linear_tolerance);
  lsParams.set("Output Frequency", 1);
  lsParams.set("Size of Krylov Subspace", 1000);

  //create linear system
  Teuchos::RCP<NOX::Epetra::Interface::Required> iReq(_interface);
  Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys;
  Teuchos::RCP<Epetra_Operator> pc;

  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
    {
      if (this->_preconditioner)
        {
          // PJNFK
          lsParams.set("Preconditioner", "User Defined");

          TrilinosPreconditioner<Number> * trilinos_pc =
            cast_ptr<TrilinosPreconditioner<Number> *>(this->_preconditioner);
          pc = Teuchos::rcp(trilinos_pc);

          Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec(_interface);
          linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iPrec, pc, x));
        }
      else
        {
          lsParams.set("Preconditioner", "None");
          //      lsParams.set("Preconditioner", "Ifpack");
          //      lsParams.set("Preconditioner", "AztecOO");

          // full jacobian
          NonlinearImplicitSystem & sys = _interface->_solver->system();
          EpetraMatrix<Number> & jacSys = *cast_ptr<EpetraMatrix<Number> *>(sys.matrix);
          Teuchos::RCP<Epetra_RowMatrix> jacMat = Teuchos::rcp(jacSys.mat());

          Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(_interface);
          linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, jacMat, x));
        }
    }
  else
    {
      // matrix free
      Teuchos::RCP<NOX::Epetra::MatrixFree> MF = Teuchos::rcp(new NOX::Epetra::MatrixFree(printParams, iReq, x));

      Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(MF);
      linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, MF, x));
    }

  //create group
  Teuchos::RCP<NOX::Epetra::Group> grpPtr = Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, x, linSys));
  NOX::Epetra::Group & grp = *(grpPtr.get());

  Teuchos::RCP<NOX::StatusTest::NormF> absresid =
    Teuchos::rcp(new NOX::StatusTest::NormF(this->absolute_residual_tolerance, NOX::StatusTest::NormF::Unscaled));
  Teuchos::RCP<NOX::StatusTest::NormF> relresid =
    Teuchos::rcp(new NOX::StatusTest::NormF(grp, this->relative_residual_tolerance));
  Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
    Teuchos::rcp(new NOX::StatusTest::MaxIters(this->max_nonlinear_iterations));
  Teuchos::RCP<NOX::StatusTest::FiniteValue> finiteval =
    Teuchos::rcp(new NOX::StatusTest::FiniteValue());
  Teuchos::RCP<NOX::StatusTest::NormUpdate> normupdate =
    Teuchos::rcp(new NOX::StatusTest::NormUpdate(this->absolute_step_tolerance));
  Teuchos::RCP<NOX::StatusTest::Combo> combo =
    Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
  combo->addStatusTest(absresid);
  combo->addStatusTest(relresid);
  combo->addStatusTest(maxiters);
  combo->addStatusTest(finiteval);
  combo->addStatusTest(normupdate);

  Teuchos::RCP<Teuchos::ParameterList> finalPars = nlParamsPtr;

  Teuchos::RCP<NOX::Solver::Generic> solver = NOX::Solver::buildSolver(grpPtr, combo, nlParamsPtr);
  NOX::StatusTest::StatusType status = solver->solve();
  this->converged = (status == NOX::StatusTest::Converged);

  const NOX::Epetra::Group & finalGroup = dynamic_cast<const NOX::Epetra::Group &>(solver->getSolutionGroup());
  const NOX::Epetra::Vector & noxFinalSln = dynamic_cast<const NOX::Epetra::Vector &>(finalGroup.getX());

  *x_epetra->vec() = noxFinalSln.getEpetraVector();
  x_in.close();

  Real residual_norm = finalGroup.getNormF();
  unsigned int total_iters = solver->getNumIterations();
  _n_linear_iterations = finalPars->sublist("Direction").sublist("Newton").sublist("Linear Solver").sublist("Output").get("Total Number of Linear Iterations", -1);

  // do not let Trilinos to deallocate what we allocated
  pc.release();
  iReq.release();

  return std::make_pair(total_iters, residual_norm);
}



template <typename T>
int
NoxNonlinearSolver<T>::get_total_linear_iterations()
{
  return _n_linear_iterations;
}



//------------------------------------------------------------------
// Explicit instantiations
template class NoxNonlinearSolver<Number>;

} // namespace libMesh



#endif // LIBMESH_TRILINOS_HAVE_NOX && LIBMESH_TRILINOS_HAVE_EPETRA