NOX Epetra Tutorial

This page describes how to write the main file for using NOX's Epetra implementation of the Group and Vector.

This example can be found in Trilinos/packages/nox/test/epetra/1Dfem/1Dfem.C

Begin by including the NOX library header files and any other headers needed by your application:

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Epetra Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_RowMatrix.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"

// User's application specific files
#include "1DfemInterface.H"
#include "1DfemPrePostOperator.H"

#include "Teuchos_ParameterList.hpp"

For convenience define the following:

using namespace std;

Begin by initializing MPI (if building a parallel version) and create the Epetra communicator for MPI. In this case autoconf defines the flag HAVE_MPI if we are building a parallel version.

int main(int argc, char *argv[])
{

  // Initialize MPI
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  // Create a communicator for Epetra objects
#ifdef HAVE_MPI
  Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
  Epetra_SerialComm Comm;
#endif

Setup some initial values.

// Get the process ID and the total number of processors
int MyPID = Comm.MyPID();
int NumProc = Comm.NumProc();

// Check verbosity level
bool verbose = false;
if (argc > 1)
  if (argv[1][0]=='-' && argv[1][1]=='v')
    verbose = true;

// Get the number of elements from the command line
int NumGlobalElements = 0;
if ((argc > 2) && (verbose))
  NumGlobalElements = atoi(argv[2]) + 1;
else if ((argc > 1) && (!verbose))
  NumGlobalElements = atoi(argv[1]) + 1;
else
  NumGlobalElements = 101;

// The number of unknowns must be at least equal to the
// number of processors.
if (NumGlobalElements < NumProc) {
  cout << "numGlobalBlocks = " << NumGlobalElements
   << " cannot be < number of processors = " << NumProc << std::endl;
  cout << "Test failed!" << std::endl;
  throw "NOX Error";
}

Create your object that derives from NOX::Epetra::Interface::Required so that nox has access to the set of nonlinear equations to solve.

// Create the interface between NOX and the application
// This object is derived from NOX::Epetra::Interface
Teuchos::RCP<Interface> interface =
  Teuchos::rcp(new Interface(NumGlobalElements, Comm));

Grab the initial guess vector and initialize it and the interface object.

// Get the vector from the Problem
Teuchos::RCP<Epetra_Vector> soln = interface->getSolution();
Teuchos::RCP<NOX::Epetra::Vector> noxSoln =
  Teuchos::rcp(new NOX::Epetra::Vector(soln,
                   NOX::Epetra::Vector::CreateView));

// Set the PDE factor (for nonlinear forcing term).  This could be specified
// via user input.
interface->setPDEfactor(1000.0);

// Set the initial guess
soln->PutScalar(1.0);

Create a parameter list and choose the options for the NOX solver. Note that the linear solver parameters in teh "Linear Solver" sublist are dependent upon what NOX::Epetra::LinearSystem object you are using. Currently NOX only comes with one concrete linear system implementation: NOX::EpetraLinearSystemAztecOO.

// Begin Nonlinear Solver ************************************

// Create the top level parameter list
Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr =
  Teuchos::rcp(new Teuchos::ParameterList);
Teuchos::ParameterList& nlParams = *(nlParamsPtr.get());

// Set the nonlinear solver method
nlParams.set("Nonlinear Solver", "Line Search Based");

// Set the printing parameters in the "Printing" sublist
Teuchos::ParameterList& printParams = nlParams.sublist("Printing");
printParams.set("MyPID", MyPID);
printParams.set("Output Precision", 3);
printParams.set("Output Processor", 0);
if (verbose)
  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 +
                           NOX::Utils::Debug +
               NOX::Utils::TestDetails +
               NOX::Utils::Error);
else
  printParams.set("Output Information", NOX::Utils::Error +
               NOX::Utils::TestDetails);

// Sublist for line search
Teuchos::ParameterList& searchParams = nlParams.sublist("Line Search");
searchParams.set("Method", "Full Step");

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

// Sublist for linear solver for the Newton method
Teuchos::ParameterList& lsParams = newtonParams.sublist("Linear Solver");
lsParams.set("Aztec Solver", "GMRES");
lsParams.set("Max Iterations", 800);
lsParams.set("Tolerance", 1e-4);
lsParams.set("Preconditioner", "New Ifpack");
lsParams.set("Preconditioner Reuse Policy", "Reuse");
lsParams.set("Max Age Of Prec", 5);

Optionally the user can define methods that will be called before and after each nonlineaer iteration and before and after each nonlinear solve. This is done by creating an object derived from NOX::Abstract::PrePostOperator.

// Add a user defined pre/post operator object
Teuchos::RCP<NOX::Abstract::PrePostOperator> ppo =
  Teuchos::rcp(new UserPrePostOperator(printing));
nlParams.sublist("Solver Options").set("User Defined Pre/Post Operator",
                   ppo);

Set the status test check option.

// Let's force all status tests to do a full check
nlParams.sublist("Solver Options").
  set("Status Test Check Type", NOX::StatusTest::Complete);

Create a Jacobian-Free Newton-Krylov method by using a NOX::Epetra::MatrixFree operator for the Jacobian.

// 2. Matrix-Free (Epetra_Operator)
Teuchos::RCP<NOX::Epetra::MatrixFree> MF =
  Teuchos::rcp(new NOX::Epetra::MatrixFree(printParams, interface,
                        noxSoln));

Create a Finite Difference operator to estimate the Jacobian matrix for preconditioning.

// 3. Finite Difference (Epetra_RowMatrix)
Teuchos::RCP<NOX::Epetra::FiniteDifference> FD =
  Teuchos::rcp(new NOX::Epetra::FiniteDifference(printParams, interface,
                          soln));

Create a linear system derived from NOX::Epetra::LinearSystem. NOX comes with a concrete implementation of the linear system class that uses AztecOO as the linear solver called NOX::Epetra::LinearSystemAztecOO. In this case we will use the constructor that specifies the Jacobian operator and a preconditioner matrix to be used with an internally constructed AztecOO preconditioner.

// Create the linear system
Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = MF;
Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec = FD;
Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys =
  Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
                            iJac, MF,
                            iPrec, FD,
                             soln));

Create the group using the linear system.

// Create the Group
NOX::Epetra::Vector initialGuess(soln, NOX::Epetra::Vector::CreateView);
Teuchos::RCP<NOX::Epetra::Group> grpPtr =
  Teuchos::rcp(new NOX::Epetra::Group(printParams,
                  iReq,
                  initialGuess,
                  linSys));
NOX::Epetra::Group& grp = *grpPtr;

Create the status tests to determine how to terminate the nonlinear solve.

// Create the convergence tests
Teuchos::RCP<NOX::StatusTest::NormF> absresid =
  Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-8));
Teuchos::RCP<NOX::StatusTest::NormF> relresid =
  Teuchos::rcp(new NOX::StatusTest::NormF(grp, 1.0e-2));
Teuchos::RCP<NOX::StatusTest::NormUpdate> update =
  Teuchos::rcp(new NOX::StatusTest::NormUpdate(1.0e-5));
Teuchos::RCP<NOX::StatusTest::NormWRMS> wrms =
  Teuchos::rcp(new NOX::StatusTest::NormWRMS(1.0e-2, 1.0e-8));
Teuchos::RCP<NOX::StatusTest::Combo> converged =
  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::AND));
converged->addStatusTest(absresid);
converged->addStatusTest(relresid);
converged->addStatusTest(wrms);
converged->addStatusTest(update);
Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
  Teuchos::rcp(new NOX::StatusTest::MaxIters(20));
Teuchos::RCP<NOX::StatusTest::FiniteValue> fv =
  Teuchos::rcp(new NOX::StatusTest::FiniteValue);
Teuchos::RCP<NOX::StatusTest::Combo> combo =
  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
combo->addStatusTest(fv);
combo->addStatusTest(converged);
combo->addStatusTest(maxiters);

Create the solver and solve the problem

// Create the solver
Teuchos::RCP<NOX::Solver::Generic> solver =
  NOX::Solver::buildSolver(grpPtr, combo, nlParamsPtr);
NOX::StatusTest::StatusType solvStatus = solver->solve();

// End Nonlinear Solver **************************************

Get the underlying solution vector in terms of an Epetra_Vector and do whatever you want with it.

// Get the Epetra_Vector with the final solution from the solver
const NOX::Epetra::Group& finalGroup =
  dynamic_cast<const NOX::Epetra::Group&>(solver->getSolutionGroup());
const Epetra_Vector& finalSolution =
  (dynamic_cast<const NOX::Epetra::Vector&>(finalGroup.getX())).
  getEpetraVector();

Exit simulation.

  return 0;
}