PETSC nonlinear solvers

Experimental html version of Parallel Programming in MPI, OpenMP, and PETSc by Victor Eijkhout. download the textbook at https:/theartofhpc.com/pcse
\[ \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: \begin{equation} \mathop{?}_x\colon f(x)=0 \end{equation} with $f\colon \bbR^n-\bbR^n$. In the special case of a linear function, \begin{equation} f(x) = Ax-b, \end{equation} 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 \begin{equation} x_{n+1} = x_n - F(x_n)\inv f(x_n) \end{equation} 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 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:

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: \begin{equation} H(x)y \approx \frac{f(x+hy)-f(x)}{h} \end{equation} 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 SNESSetJacobian  , where the Jacobian is a function of type 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 \begin{equation} u_t = G(t,u) \end{equation} call TSSetRHSFunction  . 

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

For implicit cases \begin{equation} F( t,u,u_t ) = 0 \end{equation} call TSSetIFunction 

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