PETSc - Portable, Extensible Toolkit for Scientific Computation

PETSc, pronounced PET-see (the S is silent), is a suite of data structures and routines for the scalable (parallel) solution of scientific applications modeled by partial differential equations.

Link	Web site
Language	C
License	BSD 2-clause license
Repository	GitLab

Lorenz equations

We would like to solve the Lorenz equations, a classical chaotic system given by:

\[\begin{split} \begin{cases} \dot{x} &= \sigma (y - x) \\ \dot{y} &= \rho x - y - x z \\ \dot{z} &= x y - \beta z \end{cases} \end{split}\]

with parameter \(\sigma=10\), \(\rho = 28\) and \(\beta = \frac{8}{3}\), and the initial state \((x_0, y_0, z_0) = (1,1,1)\). We solve it with a classical Runge-Kutta method of order 4, given by TSRK4 method.

To solve a problem with PETSc we first write our equation as a ODE of the form:

\[ \dot{u} = f(t, u) \]

and the user provide the function \(f\) as:

  PetscErrorCode f(TS ts, double t, Vec u, Vec du, void* params);

where u the current state of the function, t the current time and output du \(f(t,u)\), and ts the time stepper.

PetscErrorCode lorenz(TS ts, double t, Vec vec_y, Vec vec_ydot, void* params)
{
    double sigma = ((double *)params)[0];
    double rho = ((double *)params)[1];
    double beta = ((double *)params)[2];

    const double* y;
    double* ydot;

    VecGetArrayRead(vec_y, &y);
    VecGetArray(vec_ydot, &ydot);

    ydot[0] = sigma*( y[1] - y[0] );
    ydot[1] = rho*y[0] - y[1] - y[0]*y[2];
    ydot[2] = y[0]*y[1] - beta*y[2];

    VecRestoreArray(vec_ydot, &ydot);
    VecRestoreArrayRead(vec_y, &y);

    return PETSC_SUCCESS;
}

Next we need to define a time stepper and set all parameters, and give it the \(f\) function

  TS ts;
  TSCreate(PETSC_COMM_SELF, &ts);
  TSSetType(ts, TSRK);
  TSRKSetType(ts, TSRK4);

  TSSetTime(ts, 0.);
  TSSetMaxTime(ts, tf);
  TSSetTimeStep(ts, dt);

  TSSetRHSFunction(ts, NULL, lorenz, params);

now we can call the time stepper by:

  TSStep(ts);

and get current state and time with TSGetSolution TSGetTime functions. The complet time loop is

  while( tn < tf )
  {
    TSGetStepNumber(ts, &n_step);
    TSGetSolution(ts, &vec_yn);
    VecGetArray(vec_yn, &yn);

    sol[n_step][0] = tn;
    sol[n_step][1] = yn[0];
    sol[n_step][2] = yn[1];
    sol[n_step][3] = yn[2];

    VecRestoreArray(vec_yn, &yn);

    TSStep(ts);

    TSGetTime(ts, &tn);
  }

where sol is an array to store all states.

For the complet example, see lorenz.c source file.

Transport equation

In this example we would like to solve the following PDE:

\[ \partial_t u + a \partial_x u = 0 \]

with \(t>0\), on the torus \(x\in[0, 1)\), a velocity \(a=1\) and the initial condition given by a hat function:

\[\begin{split} u(0, x) = \begin{cases} x - 0.25 & \text{if } x\in[0.25, 0.5[ \\ -x + 0.75 & \text{if } x\in[0.5, 0.75[ \\ 0 & \text{else} \end{cases} \end{split}\]

We choose a first order up-wind scheme to estimate the \(x\) derivative and a forward Euler method for the time discretization given by TSRK1FE method.

We define the up-wind scheme as:

PetscErrorCode
transport( TS ts, double t, Vec vec_y, Vec vec_dy, void* params )
{
    double a  = ( (double*)params )[0];
    double dx = ( (double*)params )[1];
    int n_x   = ( (double*)params )[2];

    double const* y;
    double* dy;

    VecGetArrayRead( vec_y, &y );
    VecGetArray( vec_dy, &dy );

    dy[0] = -( fmax( a, 0. ) * ( y[0] - y[n_x - 1] ) + fmin( a, 0. ) * ( y[1] - y[0] ) ) / dx;

    for ( int i = 1; i < n_x - 1; ++i )
    {
        dy[i] = -( fmax( a, 0. ) * ( y[i] - y[i - 1] ) + fmin( a, 0. ) * ( y[i + 1] - y[i] ) ) / dx;
    }

    dy[n_x - 1] = -( fmax( a, 0. ) * ( y[n_x - 1] - y[n_x - 2] ) + fmin( a, 0. ) * ( y[0] - y[n_x - 1] ) ) / dx;

    VecRestoreArray( vec_dy, &dy );
    VecRestoreArrayRead( vec_y, &y );

    return PETSC_SUCCESS;
}

The time loop is the same as for Lorenz equation.

For the complet example, see transport.c source file.

Arenstorf orbit

The Arenstorf orbit problem is a classical problem to test adaptive time step methods:

\[\begin{split} \begin{cases} \ddot{x} &= x + 2\dot{y} - \frac{1-\mu}{r_1^3}(x+\mu) - \frac{\mu}{r_2^3}(x-1+\mu) \\ \ddot{y} &= y - 2\dot{x} - \frac{}{1-\mu}{r_1^3}y - \frac{\mu}{r_2^3}y \end{cases} \end{split}\]

with initial condition \((x,\dot{x},y,\dot{y})=(0.994, 0, 0, -2.001585106)\), \(r_1\) and \(r_2\) given by

\[ r_1 = \sqrt{(x+\mu)^2 + y^2},\quad r_2 = \sqrt{(x-1+\mu)^2 + y^2} \]

and with parameter \(\mu = 0.012277471\).

First of all, we need to rewrite this problem into a first order derivative equation in time

\[\begin{split} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} = \begin{pmatrix} x \\ y \\ \dot{x} \\ \dot{y} \end{pmatrix}, \qquad \begin{cases} \dot{y}_1 = y_3 \\ \dot{y}_2 = y_4 \\ \dot{y}_3 = y_1 + 2y_4 - \frac{1-\mu}{r_1^3}(y_1 + \mu) - \frac{\mu}{r_2^3}(y_1-1+\mu) \\ \dot{y}_4 = y_2 - 2y_3 - \frac{1-\mu}{r_1^3}y_2 - \frac{\mu}{r_2^3}y_2 \\ \end{cases} \end{split}\]

We define this system as:

PetscErrorCode arenstorf(TS ts, double t, Vec vec_y, Vec vec_ydot, void* params)
{
    double mu = *(double *)params;

    const double* y;
    double* ydot;

    VecGetArrayRead(vec_y, &y);
    VecGetArray(vec_ydot, &ydot);

    double const y1 = y[0];
    double const y2 = y[1];
    double const y3 = y[2];
    double const y4 = y[3];

    double const r1 = sqrt( ( y1 + mu ) * ( y1 + mu ) + y2 * y2 );
    double const r2 = sqrt( ( y1 - 1. + mu ) * ( y1 - 1. + mu ) + y2 * y2 );

    ydot[0] = y3;
    ydot[1] = y4;
    ydot[2] = y1 + 2. * y4 - ( 1. - mu ) * ( y1 + mu ) / ( r1 * r1 * r1 ) - mu * ( y1 - 1. + mu ) / ( r2 * r2 * r2 );
    ydot[3] = y2 - 2. * y3 - ( 1. - mu ) * y2 / ( r1 * r1 * r1 ) - mu * y2 / ( r2 * r2 * r2 );

    VecRestoreArray(vec_ydot, &ydot);
    VecRestoreArrayRead(vec_y, &y);

    return PETSC_SUCCESS;
}

We solve this example with given method TSRK5DP which is the method RK8(7) 13M in [PD81] (mainly call DOPRI8). For this we need to set the integrator and settings for adaptive time step method

  TS ts;
  TSCreate(PETSC_COMM_SELF, &ts);
  TSSetType(ts, TSRK);
  TSRKSetType(ts, TSRK5DP);
  TSSetTolerances(ts, 1e-5, NULL, 1e-5, NULL);

  TSSetTime(ts, 0.);
  TSSetMaxTime(ts, tf);
  TSSetTimeStep(ts, dt);

  TSSetRHSFunction(ts, NULL, arenstorf, &mu);

then the time loop becomes

  while( tn < tf )
  {
    TSGetStepNumber(ts, &n_step);
    TSGetTimeStep(ts, &dt);
    TSGetSolution(ts, &vec_yn);
    VecGetArray(vec_yn, &yn);

    sol[n_step][0] = tn;
    sol[n_step][1] = yn[0];
    sol[n_step][2] = yn[1];
    sol[n_step][3] = yn[2];
    sol[n_step][4] = yn[3];
    sol[n_step][5] = dt;

    VecRestoreArray(vec_yn, &yn);

    TSStep(ts);

    TSGetTime(ts, &tn);
  }

For the complet example, see arenstorf.c source file.