xx Solving Nonlinear Elliptic Equations: Calling PETSc from Cactus

Solving Nonlinear Elliptic Equations: Calling PETSc from Cactus

Erik Schnetter schnetter@uni-tuebingen.de

August 28, 2001

Abstract

The Cactus thorn TATPETSc provides a simple interface to the SNES interface of PETSc. SNES (“Simple Nonlinear Elliptic Solver”) is a set of efficient parallel solvers for sparse matrices, i.e.xfor discretised problems with small stencils. The main task of TATPETSc is in handling the mismatch in the different parallelisation models.

Contents

1 Introduction 2 Nonlinear Elliptic Equations 3 Solving a Nonlinear Elliptic Equation 4 Initial data 5 Boundary Conditions 6 The Solution, Or Not 7 Excision 8 Common PETSc options 9 Installing PETSc 10 Parameters 11 Interfaces 12 Schedule

1 Introduction PETSc (“The Portable, Extensible Toolkit for Scientific Computation”) [1] is a library that, among other things, solves nonlinear elliptic equations. It is highly configurable, and using it is a nontrivial exercise. It is therefore convenient to create wrappers for this library that handle the more common cases. Although this introduces a layer that is a black box to the PETSc novice, it at the same time enables this PETSc novice to use this library, which would otherwise not be possible.

At the moment there exist three wrapper routines. Two of them handle initialisation and shutting down the library at programme startup and shutdown time. The third solves a nonlinear elliptic equation.

2 Nonlinear Elliptic Equations Nonlinear elliptic equations can be written in the form where is an -dimensional vector, and is function with an -dimensional vector as result. The vector is the unknown to solve for. The function has to be given together with initial data for the unknown .

In the case of discretised field equations, the vector contains all grid points. The function then includes the boundary conditions.

Consider as an example the Laplace equation in three dimensions on a grid with points. In that case, the unknown is , with the dimension . The function is given by . Remember that, once a problem is discretised, derivative operators are algebraic expressions. That means that calculating does not involve taking analytic (“real”) derivatives, but can rather be written as linear function (matrix multiplication).

3 Solving a Nonlinear Elliptic Equation The wrapper routine that solves a nonlinear elliptic equation is

void TATPETScxsolve (
                                            const cGH *cctkGH,
                                               const int *var,
                                               const int *val,
                                                    int nvars,
                                            int optionsxtable,
        void (*fun) (cGH *cctkGH, int optionsxtable, void *userdata),
        void (*bnd) (cGH *cctkGH, int optionsxtable, void *userdata),

                                             void *data);

There is currently no Fortran wrapper for this routine. This routine takes the following arguments:

cctkGH

a pointer to the Cactus grid hierarchy

var

the list of indices of the grid variables to solve for

val

the list of indices of the grid variables that contain the function value

nvars

the number of variables to solve for, which is the same as the number of function values

optionsxtable

a table with additional options (see below)

fun

the routine that calculates the function values from the variables

bnd

the routine that applies the boundary conditions to the variables

data

data that are passed through unchanged to the callback routines

The options table can have the following elements:

CCTKxINT periodic[dim]

flags indicating whether the grid is periodic in the corresponding direction. ( is the number of dimensions of the grid variables.) The default values are 0.

CCTKxINT solvebnds[2*dim]

flags indicating whether the grid points on the corresponding outer boundaries should be solved for (usually not). The default values are 0.

CCTKxINT stencilxwidth

the maximum stencil size used while calculating the function values from the variables (should be as small as possible). The default value is 1.

CCTKxFNxPOINTER jacobian
The real type of jacobian must be void (*jacobian)(cGH *cctkGH,

void *data). This is either a routine that calculates the Jacobian directly, or 0 (zero). The default is 0.

CCTKxFNxPOINTER getxcoloring
The real type of getxcoloring must be void (*getxcoloring)(DA
da, ISColoring *iscoloring, void *userdata). This is either a

routine that calculates the colouring for calculating the Jacobian, or 0 (zero). You need to pass a nonzero value only if you have very special boundary conditions. The default is 0.

A few remarks:

In order to be able to call this routine directly, it is necessary to inherit from TATPETSc in the thorn where this routine is called.

Increasing stencilxwidth from 1 to 2 increases the run time by about a factor of 5. In general you want to be able to set stencilxwidth to 1. Note that calculating does normally not require upwind derivatives with their larger size. stencilxwidth does not have to be equal to the number of ghost zones.

One sure way to speed up solving nonlinear elliptic equations is to provide an explicit function that calculates the Jacobian. While such a function is in principle straightforward to write, it is very (very!) tedious to do so. If you pass 0 for this function, the Jacobian is evaluated numerically. This is about one order of magnitude slower, but is also a lot less work.

4 Initial data A linear elliptic equation does not need initial data. However, a nonlinear elliptic equation does. It may have several solutions, and the initial data select between these solutions. The initial data have to be put into the variables before the solver is called. The function values can remain undefined.

5 Boundary Conditions There are two functions that you have to apply boundary conditions to. You need to apply boundary conditions to the right hand side, i.e.x, and to the solution, i.e.x.

The routine bnd has to apply boundary conditions to at exactly those boundaries where solvebnds is false. (The boundaries where solvebnds is true have already been determined by the solver.) This includes imposing the symmetry boundary conditions. Not applying a boundary condition on the outer boundaries is equivalent to a Dirichlet boundary condition, with the boundary values given in the initial data .

After calculating , you need to apply the necessary boundary conditions to as well. You have to do this in the routine fun that calculates .

6 The Solution, Or Not When the solver returns, the result is available in . contains the corresponding function value, which should be close to zero.

It is possible that the solver doesn’t converge. In this case, will not be close to zero.

It is possible (by the way of programming error) to make certain grid points in independent of , or to ignore certain grid points in while calculating . This leads to a singular matrix when the nonlinear solver calls a linear solver as a substep. Such a system is ill-posed and cannot be solved.

7 Excision An excision boundary leads to a certain number of grid points that should not be solved for. In order to avoid a singular matrix, it is still necessary to impose a condition on these grid points. Assuming that you want a Dirichlet-like condition for these grid points, I suggest where are the initial data that you have to save someplace. Note that you have to impose this (or a different) condition not only onto the boundary points, but onto all excised points, i.e.xall points that are not solved for.

(The above condition satisfies for , which will be the solution for these grid points. Setting does not work because it leads to a singular matrix, as outlined above.)

8 Common PETSc options Options Database Keys (from the PETSc documentation)

-snesxtype type

ls, tr, umls, umtr, test

-snesxstol

convergence tolerance in terms of the norm of the change in the solution between steps

-snesxatol atol

absolute tolerance of residual norm

-snesxrtol rtol

relative decrease in tolerance norm from initial

-snesxmaxxit maxxit

maximum number of iterations

-snesxmaxxfuncs maxxfuncs

maximum number of function evaluations

-snesxtrtol trtol

trust region tolerance

-snesxnoxconvergencextest

skip convergence test in nonlinear or minimization solver; hence iterations will continue until maxxit or some other criterion is reached. Saves expense of convergence test

-snesxmonitor

prints residual norm at each iteration

-snesxvecmonitor

plots solution at each iteration

-snesxvecmonitorxupdate

plots update to solution at each iteration

-snesxxmonitor

plots residual norm at each iteration

-snesxfd

use finite differences to compute Jacobian; very slow, only for testing

-snesxmfxkspxmonitor

if using matrix-free multiply then print at each KSP iteration

9 Installing PETSc Before you can use TATPETSc, you have to install PETSc. PETSc comes with extensive documentation for that.

In order to be able to use PETSc with Cactus, you have to give certain options when you configure your Cactus applications. To do so, create an options file containing the following options.

SYSxINCxDIRS /usr/include/petsc

LIBDIRS /usr/X11R6/lib

LIBS crypt petscfortran petscts petscsnes petscsles petscdm petscmat petscvec petsc lapack blas mpe mpich X11 g2c z

Replace /usr/include/petsc with the corresponding directory on your machine.

If you have other packages that need other options, then you have to combine these options manually. Even if the other options would normally be selected automatically, selecting the PETSc options manually will override the other options.

Do not forget to activate MPI. PETSc needs MPI to run.

References

  [1]xxxPETSc: http://www-fp.mcs.anl.gov/petsc/

10 Parameters

11 options 10 Scope: privatex01 STRING

11 Description: Command line options for PETSc

11 Range 10 x01 Default: (none)
11 .*  10 no restriction

11

11 verbose  10 Scope: privatex01  BOOLEAN

11 Description: Produce log output while running

11            10  x01  Default: no

11

11 veryverbose 10 Scope: privatex01  BOOLEAN

11 Description: Produce much log output while running

11            10  x01  Default: no

11

11 Interfaces

General Implements:

tatpetsc

Inherits:

tatelliptic

Uses header:

TATelliptic.h

12 Schedule This section lists all the variables which are assigned storage by thorn CactusElliptic/TATPETSc. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage NONE

Scheduled Functions

CCTKxSTARTUP x tatpetscxinitialize x initialise petsc

11 x 11 After:  11 driverxstartup
11 x 11 Language: 11 c
11 x 11 Type:  11 function

11

CCTKxSHUTDOWN x tatpetscxfinalize x finalise petsc x

11 x 11 Language: 11 c
11 x 11 Type:  11 function