PyTrilinos.NOX

Python versions of C++ methods that expect Teuchos::ParameterList arguments accept either Teuchos.ParameterList objects or python dictionaries.

Since python objects are already reference counted, the Teuchos::RCP used in the C++ version of NOX is hidden from python users. Python versions of C++ methods that accept Teuchos::RCP<object> arguments accept raw object arguments in python.

The following NOX namespaces are suported in python:

NOX
NOX.Abstract
NOX.Epetra
NOX.Epetra.Interface
NOX.Solver
NOX.StatusTest

The only NOX Interface class that has a python wrapper is the Epetra interface. As with the C++ version of NOX, you define a nonlinear problem by declaring a python class that inherits from NOX.Epetra.Interface.Required. Optionally, your class may also inherit from NOX.Epetra.Interface.Jacobian and/or NOX.Epetra.Interface.Preconditioner.

Your constructor must call the constructors of its NOX.Epetra.Interface base classes:

class MyProblem(NOX.Epetra.Interface.Required):
  def __init__(self):
    NOX.Epetra.Interface.Required.__init__(self)

This is necessary because the underlying callback mechanism needs to be properly initialized. The nonlinear function for your class is implemented by defining a computeF() method:

def computeF(self, x, F, flag):
  "Required implementation of computeF() method"
  ...

Arguments x and F are passed from the underlying C++ code as Epetra_Vector objects and then converted to numpy-hybrid Epetra.Vector objects when they are passed into your python computeF() function. Therefore, you can use the numpy.ndarray interface on these arguments, such as the shape attribute, or slice indexing.

The flag argument is an integer that NOX uses to inform computeF() why it is being called. You can use this flag to alter how computeF() behaves. See the C++ NOX documentation for details.

We are defining a somewhat complicated python-to-C++-to-python callback scheme here. Unfortunately, error-handling in such an environment is not as rubust as in a single-language environment. Specifically, if your computeF() (accidentally) raises a python exception, all you may see is:

terminate called after throwing an instance of 'Swig::DirectorMethodException'
Abort trap

For this reason, it is a good idea to wrap your computations in a try block:

class MyProblem(NOX.Epetra.Interface.Required):

  def __init__(self):
    NOX.Epetra.Interface.Required.__init__(self)

  def computeF(self, x, F, flag):
    "Required implementation of computeF() method"
    try:
      F[:] = nonlinear_func(x)  # Whatever this may be...
    except Exception, e:
      print "Python exception raised in MyProblem.computeF:"
      print e
      return False
    return True

This will print the python exception if one is raised, and tells NOX that the computation failed. Remember to return a boolean indicating success or failure of the computeF() method.

NOX needs to be able to compute the result of the Jacobian of F(x) on a given vector from the same space as x. The NOX.Epetra.Interface.Required class and computeF() method are sufficient to estimate this product if you are willing to use the NOX.Epetra.FiniteDifference class. If you want to use an algebraic preconditioner based on an approximation to the Jacobian, you will want to use map coloring to greatly improve the efficiency, specifically the NOX.Epetra.FiniteDifferenceColoring class. See exNOX_1Dfem.py in the example directory of the PyTrilinos package for an example of this type of implementation.

You solve your problem by first creating an instance of your MyProblem class and an Epetra.Vector solution (using some appropriate constructor):

problem = MyProblem()
soln    = Epetra.Vector(...)

Since we do not inherit from Jacobian or Preconditioner, we need to define a matrix-free linear system for NOX to use:

mf      = NOX.Epetra.MatrixFree({ }, problem, soln)
fdc     = NOX.Epetra.FiniteDifferenceColoring({ }, { }, soln,
                                              problem.getGraph())
linSys  = NOX.Epetra.LinearSystemAztecOO({ }, { }, mf, mf, fdc, fdc,
                                         soln)

Note that these calls assume that your MyProblem class defines a getGraph() method that returns an Epetra.CrsGraph that defines the sparsity/coupling of the nonlinear problem. Also, the empty python dictionaries result in default parameter specifications. Now we can create a NOX.Epetra.Group and a NOX.Solver.Manager:

group   = NOX.Epetra.Group({ }, problem, soln, linSys)
solver  = NOX.Solver.Manager(group, statusTest, { })

Here we assume statusTest is a properly constructed NOX.StatusTest object. Now that everything is specified, we can solve the problem with:

status  = solver.solve()

Simple, no?