PETSC nonlinear solvers

Experimental html version of downloadable textbook, see http://www.tacc.utexas.edu/~eijkhout/istc/istc.html
\[ \newcommand\inv{^{-1}}\newcommand\invt{^{-t}} \newcommand\bbP{\mathbb{P}} \newcommand\bbR{\mathbb{R}} \newcommand\defined{ \mathrel{\lower 5pt \hbox{${\equiv\atop\mathrm{\scriptstyle D}}$}}} \] 36.1 : Nonlinear systems
36.1.1 : Basic setup
36.2 : Time-stepping
Back to Table of Contents

36 PETSC nonlinear solvers

36.1 Nonlinear systems

crumb trail: > petsc-nonlinear > Nonlinear systems

Nonlinear system solving means finding the zero of a general nonlinear function, that is: \[ \mathop{?}_x\colon f(x)=0 \] with $f\colon \bbR^n-\bbR^n$. In the special case of a linear function, \[ f(x) = Ax-b, \] we solve this by any of the methods in chapter~ PETSc solvers .

The general case can be solved by a number of methods, foremost Newton's method , which iterates \[ x_{n+1} = x_n - F(x_n)\inv f(x_n) \] where $F$ is the Hessian $F_{ij}=\partial f_i/\partial x_j$.

You see that you need to specify two functions that are dependent on your specific problem: the objective function itself, and its Hessian.

36.1.1 Basic setup

crumb trail: > petsc-nonlinear > Nonlinear systems > Basic setup

The PETSc nonlinear solver object is of type SNES : `simple nonlinear equation solver'. As with linear solvers, we create this solver on a communicator, set its type, incorporate options, and call the solution routine

#include "petscsnes.h"
PetscErrorCode  SNESSolve(SNES snes,Vec b,Vec x)

Collective on SNES

Input Parameters
snes - the SNES context
b    - the constant part of the equation F(x) = b, or NULL to use zero.
x    - the solution vector.

SNESSolve :

Vec value_vector,solution_vector;
/* vector creation code missing */
SNES solver;
SNESCreate( comm,&solver );
SNESSetFunction( solver,value_vector,formfunction, NULL );
SNESSetFromOptions( solver );
SNESSolve( solver,NULL,solution_vector );

The function has the type

PetscErrorCode formfunction(SNES,Vec,Vec,void*)

where the parameters are:

  • the solver object, so that you can access to its internal parameters
  • the $x$ value at which to evaluate the function
  • the result vector $f(x)$ for the given input
  • a context pointer for further application-specific information.

Example:

PetscErrorCode evaluation_function( SNES solver,Vec x,Vec fx, void *ctx ) {
  const PetscReal *x_array;
  PetscReal *fx_array;
  VecGetArrayRead(fx,&fx_array);
  VecGetArray(x,&x_array);
  for (int i=0; i<localsize; i++)
    fx_array[i] = pointfunction( x_array[i] );
  VecRestoreArrayRead(fx,&fx_array);
  VecRestoreArray(x,&x_array);
};

Comparing the above to the introductory description you see that the Hessian is not specified here. An analytic Hessian can be dispensed with if you instruct PETSc to approximate it by finite differences: \[ H(x)y \approx \frac{f(x+hy)-f(x)}{h} \] with $h$ some finite diference. The commandline option snes_fd forces the use of this finite difference approximation. However, it may lead to a large number of function evaluations. The option snes_fd_color applies a coloring to the variables, leading to a drastic reduction in the number of function evaluations.

If you can form the analytic Jacobian / Hessian, you can specify it with

#include "petscsnes.h"
PetscErrorCode  SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx)

Logically Collective on SNES

Input Parameters
snes - the SNES context
Amat - the matrix that defines the (approximate) Jacobian
Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
J    - Jacobian evaluation routine (if NULL then SNES retains any previously set value)
ctx  - [optional] user-defined context for private data for the Jacobian evaluation routine
SNESSetJacobian , where the Jacobian is a function of type
#include "petscsnes.h"
PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx);

Collective on snes

Input Parameters
x   - input vector, the Jacobian is to be computed at this value
ctx - [optional] user-defined Jacobian context

Output Parameters
Amat - the matrix that defines the (approximate) Jacobian
Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
SNESJacobianFunction .

Specifying the Jacobian:

Mat J;
ierr = MatCreate(comm,&J); CHKERRQ(ierr);
ierr = MatSetType(J,MATSEQDENSE); CHKERRQ(ierr);
ierr = MatSetSizes(J,n,n,N,N); CHKERRQ(ierr);
ierr = MatSetUp(J); CHKERRQ(ierr);
ierr = SNESSetJacobian(solver,J,J,&Jacobian,NULL); CHKERRQ(ierr);

36.2 Time-stepping

crumb trail: > petsc-nonlinear > Time-stepping

For cases \[ u_t = G(t,u) \] call TSSetRHSFunction .

#include "petscts.h"
PetscErrorCode TSSetRHSFunction
   (TS ts,Vec r,
    PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),
    void *ctx);

For implicit cases \[ F( t,u,u_t ) = 0 \] call TSSetIFunction

#include "petscts.h"
PetscErrorCode TSSetIFunction
   (TS ts,Vec r,TSIFunction f,void *ctx)
Back to Table of Contents