PFASST¶

The Parallel Full Approximation Scheme in Space and Time (PFASST) algorithm is a method for parallelizing ODEs and PDEs in time. The maths behind the PFASST algorithm are breifly described on the Maths page.

This work was supported by the Director, DOE Office of Science, Office of Advanced Scientific Computing Research, Office of Mathematics, Information, and Computational Sciences, Applied Mathematical Sciences Program, under contract DE-SC0004011. This work is currently authored by Michael L. Minion and Matthew Emmett. Contributions are welcome – please contact Matthew Emmett.

PyPFASST¶

PyPFASST is a Python implementation of the PFASST algorithm.

Main parts of the documentation

Download - download and installation instructions.
Tutorial - getting started and basic usage.
Overview - design and interface overview.
Reference - reference and API documentation.

Maths¶

The PFASST algorithm is, in some ways, similar to the parareal method for the temporal parallelization of ODEs and PDEs. These methods can be roughly described in terms of two numerical approximation methods, denoted here by \(\mathcal{G}\) and \(\mathcal{F}\). Both \(\mathcal{G}\) and \(\mathcal{F}\) propagate an initial value \(U_n \approx u(t_n)\) by approximating the solution to

\[u(t) = u_0 + \int_0^t f(\tau,u(\tau)) \,d\tau\]

from \(t_n\) to \(t_{n+1}\) (that is, \(\dot{u} = f(u,t)\)). For the methods to be efficient, it must be the case that the \(\mathcal{G}\) propagator is computationally less expensive than the \(\mathcal{F}\) propagator, and hence \(\mathcal{G}\) is typically a low-order method. Since the overall accuracy of the methods is limited by the accuracy of the \(\mathcal{F}\) propagator, \(\mathcal{F}\) is typically higher-order and in addition may use a smaller time step than \(\mathcal{G}\). For these reasons, \(\mathcal{G}\) is referred to as the coarse propagator and \(\mathcal{F}\) the fine propagator.

For PDEs, it may be the case that the coarse and fine propagators differ only in their level of refinement (in the spatial and/or temporal directions).

The parareal method begins by computing a first approximation \(U_{n+1}^0\) for \(n = 0 \ldots N-1\) where \(N\) is the number of time steps, often performed with the coarse propagator:

\[U_{n+1}^0 = \mathcal{G}(t_{n+1}, t_{n}, U_n^0).\]

The parareal method then proceeds iteratively, alternating between the parallel computation of \(\mathcal{F}(t_{n+1},t_n,U_n^k)\) and an update of the initial conditions at each processor of the form

\[U_{n+1}^{k+1} = \mathcal{G}(t_{n+1}, t_n, U_n^{k+1}) + \mathcal{F}(t_{n+1}, t_n, U_n^k) - \mathcal{G}(t_{n+1}, t_n, U_n^{k})\]

for \(n = 0 \ldots N-1\). That is, the fine propagator is used to refine the solution in each time slice in parallel, while the coarse propagator is used to propagate the refinements performed by the fine propagator through time to later processors.

The PFASST method operaters in a similar manner to the parareal method, but it employs the SDC time-integration technique in a novel way to improve the parallel efficiency of the method. However, the general structure and roles of the fine and coarse propagators are the same as the parareal method.

Furthermore, the PFASST algorithm adds corrections to the coarse propagator in a manner reminiscent of FAS corrections for non-linear multi-grid methods.

References¶

Download¶

The source code for PyPFASST is hosted on GitHub. Please see the PyPFASST project page to browse the code, submit bug reports, contribute etc.

To download the latest development version of PyPFASST using git:

$ git clone git://github.com/memmett/PyPFASST.git

If you have already done this and just need to grab the latest version:

$ cd PyPFASST
$ git pull

If you don’t have git, you can download the latest version of PyPFASST as a

tarball: https://github.com/memmett/PyPFASST/tarball/master
zipball: https://github.com/memmett/PyPFASST/zipball/master.

Installation¶

PyPFASST requires two Python packages: mpi4py and NumPy; and an MPI library. To install these packages on Debian or Ubuntu:

$ sudo apt-get install python-dev python-numpy mpich2 libmpich2-dev

Next, you need to install mpi4py. The easiest way to do this is with PIP:

$ sudo apt-get install python-pip
$ sudo pip install mpi4py

Tutorial¶

Advection/Diffusion¶

To run the advection/diffustion (AD) example:

$ cd PyPFASST
$ export PYTHONPATH=$PYTHONPATH:$PWD
$ cd examples/advection
$ mpirun -n 4 python main.py

This solves the 1d AD equation using a V cycle with 3 PFASST levels (5, 3, and 2 nodes) and 4 processors. The logarithm of the maximum absolute error (compared to an exact solution) is echoed after each SDC sweep.

Please skim over the overview documentation to get a jist of how PyPFASST is used, then consider the main.py script of the AD example:

"""Solve the advection/diffusion equation with PyPFASST."""

# Copyright (c) 2011, Matthew Emmett.  All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
#   1. Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
#
#   2. Redistributions in binary form must reproduce the above
#      copyright notice, this list of conditions and the following
#      disclaimer in the documentation and/or other materials provided
#      with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.


from mpi4py import MPI

import argparse

import pfasst
import pfasst.imex

from ad import *


######################################################################
# options

parser = argparse.ArgumentParser(
    description='solve the advection/diffusion equation')
parser.add_argument('-d',
                    type=int,
                    dest='dim',
                    default=1,
                    help='number of dimensions, defaults to 1')
parser.add_argument('-n',
                    type=int,
                    dest='steps',
                    default=MPI.COMM_WORLD.size,
                    help='number of time steps, defaults to number of mpi processes')
parser.add_argument('-l',
                    type=int,
                    dest='nlevs',
                    default=3,
                    help='number of levels, defaults to 3')
options = parser.parse_args()


###############################################################################
# config

comm  = MPI.COMM_WORLD
nproc = comm.size

dt   = 0.01
tend = dt*options.steps

N = 1024
D = options.dim

nnodes = [ 9, 5, 3 ]


###############################################################################
# init pfasst

pf = pfasst.PFASST()
pf.simple_communicators(ntime=nproc, comm=comm)

for l in range(options.nlevs):
  F = AD(shape=D*(N,), refinement=2**l, dim=D)
  SDC = pfasst.imex.IMEXSDC('GL', nnodes[l])
  pf.add_level(F, SDC, interpolate, restrict)

if len(pf.levels) > 1:
  pf.levels[-1].sweeps = 2


###############################################################################
# add hooks

def echo_error(level, state, **kwargs):
  """Compute and print error based on exact solution."""

  if level.feval.burgers:
    return

  y1 = np.zeros(level.feval.shape)
  level.feval.exact(state.t0+state.dt, y1)

  err = np.log10(abs(level.qend-y1).max())

  print 'step: %03d, iteration: %03d, position: %d, level: %02d, error: %f' % (
    state.step, state.iteration, state.cycle, level.level, err)

pf.add_hook(0, 'post-sweep', echo_error)


###############################################################################
# create initial condition and run

F  = AD(shape=D*(N,), dim=D)
q0 = np.zeros(F.shape)
F.exact(0.0, q0)

pf.run(q0=q0, dt=dt, tend=tend)

On line 38 we import the AD class and the spatial interpolate and restrict functions for this example from the ad.py file:

"""Solve various advection/diffusion type equations with PyPFASST."""

# Copyright (c) 2011, Matthew Emmett.  All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
#   1. Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
#
#   2. Redistributions in binary form must reproduce the above
#      copyright notice, this list of conditions and the following
#      disclaimer in the documentation and/or other materials provided
#      with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.


import math

import numpy as np
import numpy.fft as fft

import pfasst.imex


###############################################################################
# define AD level

class AD(pfasst.imex.IMEXFEval):
  """IMEX FEval class for the adv/diff equation (or viscous Burgers)."""

  def __init__(self, shape=None, refinement=1,
               dim=1, L=1.0, nu=0.005, v=1.0, t0=1.0,
               burgers=False, verbose=False, **kwargs):

    super(AD, self).__init__()

    N = shape[0] / refinement

    self.shape = dim*(N,)
    self.size  = N**dim
    self.N     = N
    self.L     = L
    self.v     = v
    self.nu    = nu
    self.t0    = t0
    self.dim   = dim
    self.burgers = burgers


    # frequencies = 2*pi/L * (wave numbers)
    K = 2*math.pi/L * fft.fftfreq(N) * N

    if verbose:
      print 'building operators...'

    # operators
    if dim == 1:

      sgradient = K*1j
      laplacian = -K**2

    elif dim == 2:

      sgradient = np.zeros(self.shape, dtype=np.complex128)
      laplacian = np.zeros(self.shape, dtype=np.complex128)

      for i in xrange(N):
        for j in xrange(N):
          laplacian[i,j] = -(K[i]**2 + K[j]**2)
          sgradient[i,j] = K[i] * 1j + K[j] * 1j

    elif dim == 3:

      sgradient = np.zeros(self.shape, dtype=np.complex128)
      laplacian = np.zeros(self.shape, dtype=np.complex128)

      for i in range(N):
        for j in range(N):
          for k in range(N):
            laplacian[i,j,k] = -(
              K[i]**2 + K[j]**2 + K[k]**2)
            sgradient[i,j,k] = (
              K[i] * 1j + K[j] * 1j + K[k] * 1j)

    else:
      raise ValueError, 'dimension must be 1, 2, or 3'

    self.sgradient  = sgradient
    self.laplacian  = laplacian

    # spectral interpolation masks
    if verbose:
      print 'building interpolation masks...'

    self.full, self.half = spectral_periodic_masks(dim, N)


  def f1_evaluate(self, u, t, f1, **kwargs):
    """Evaluate explicit piece."""

    z = fft.fftn(u)

    z_sgrad = self.sgradient * z
    u_sgrad = np.real(fft.ifftn(z_sgrad))

    if self.burgers:
      f1[...] = - (u * u_sgrad)
    else:
      f1[...] = - (self.v * u_sgrad)


  def f2_evaluate(self, y, t, f2, **kwargs):
    """Evaluate implicit piece."""

    z = fft.fftn(y)
    z = self.nu * self.laplacian * z
    u = np.real(fft.ifftn(z))

    f2[...] = u


  def f2_solve(self, rhs, y, t, dt, f2, **kwargs):
    """Solve and evaluate implicit piece."""

    # solve (rebuild operator every time, as dt may change)
    invop = 1.0 / (1.0 - self.nu*dt*self.laplacian)

    z = fft.fftn(rhs)
    z = invop * z

    y[...] = np.real(fft.ifftn(z))

    # evaluate
    z = self.nu * self.laplacian * z

    f2[...] = np.real(fft.ifftn(z))


  def exact(self, t, q):
    """Exact solution (for adv/diff equation)."""

    if self.burgers:
      raise ValueError, 'exact solution not known for non-linear case'

    dim  = self.dim
    L    = self.L
    v    = self.v
    nu   = self.nu
    t0   = self.t0

    q1 = np.zeros(self.shape, q.dtype)

    images = range(-1,2)
    #images = [0]

    if dim == 1:

      nx, = self.shape

      for ii in images:
        for i in range(nx):
          x = L*(i-nx/2)/nx + ii*L - v*t
          q1[i] += ( (4.0*math.pi*nu*(t+t0))**(-0.5)
                     * math.exp(-x**2/(4.0*nu*(t+t0))) )

    elif dim == 2:

      nx, ny = self.shape

      for ii in images:
        for jj in images:

          for i in range(nx):
            x = L*(i-nx/2)/nx + ii*L - v*t
            for j in range(ny):
              y = L*(j-ny/2)/ny + jj*L - v*t

              q1[i,j] += ( (4.0*math.pi*nu*(t+t0))**(-1.0)
                           * math.exp(-(x**2+y**2)/(4.0*nu*(t+t0))) )


    elif dim == 3:

      nx, ny, nz = self.shape

      for ii in images:
        for jj in images:
          for kk in images:

            for i in range(nx):
              x = L*(i-nx/2)/nx + ii*L - v*t
              for j in range(ny):
                y = L*(j-ny/2)/ny + jj*L - v*t
                for k in range(nz):
                  z = L*(j-nz/2)/nz + kk*L - v*t

                  q1[i,j,k] += ( (4.0*math.pi*nu*(t+t0))**(-1.5)
                                 * math.exp(-(x**2+y**2+z**2)/(4.0*nu*(t+t0))) )


    q[...] = q1


###############################################################################
# define interpolator and restrictor

def interpolate(yF, yG, fevalF=None, fevalG=None,
                dim=1, xrat=2, interpolation_order=-1, **kwargs):
  """Interpolate yG to yF."""

  if interpolation_order == -1:

    zG = fft.fftn(yG)
    zF = np.zeros(fevalF.shape, zG.dtype)

    zF[fevalF.half] = zG[fevalG.full]

    yF[...] = np.real(2**dim*fft.ifftn(zF))

  elif interpolation_order == 2:

    if dim != 1:
      raise NotImplementedError

    yF[0::xrat] = yG
    yF[1::xrat] = (yG + np.roll(yG, -1)) / 2.0

  elif interpolation_order == 4:

    if dim != 1:
      raise NotImplementedError

    yF[0::xrat] = yG
    yF[1::xrat] = ( - np.roll(yG,1)
                    + 9.0*yG
                    + 9.0*np.roll(yG,-1)
                    - np.roll(yG,-2) ) / 16.0

  else:
    raise ValueError, 'interpolation order must be -1, 2 or 4'


def restrict(yF, yG, fevalF=None, fevalG=None,
             dim=1, xrat=2, **kwargs):
  """Restrict yF to yG."""

  if yF.shape == yG.shape:
    yG[:] = yF[:]

  elif dim == 1:
    yG[:] = yF[::xrat]

  elif dim == 2:
    y = np.reshape(yF, fevalF.shape)
    yG[...] = y[::xrat,::xrat]

  elif dim == 3:
    y = np.reshape(yF, fevalF.shape)
    yG[...] = y[::xrat,::xrat,::xrat]


###############################################################################
# helpers

from itertools import product

def spectral_periodic_masks(dim, N, **kwargs):
  """Spectral interpolation masks."""

  if dim == 1:
    mask_half = np.zeros(N, dtype=np.bool)
    mask_full = np.zeros(N, dtype=np.bool)

    full = range(N)
    full.remove(N/2)
    mask_full[full] = True

    half = range(N/4) + range(3*N/4+1, N)
    mask_half[half] = True

  elif dim == 2:
    mask_half = np.zeros((N, N), dtype=np.bool)
    mask_full = np.zeros((N, N), dtype=np.bool)

    full = range(N)
    full.remove(N/2)
    for i in product(full, full):
      mask_full[i] = True

    half = range(N/4) + range(3*N/4+1, N)
    for i in product(half, half):
      mask_half[i] = True

  elif dim == 3:
    mask_half = np.zeros((N, N, N), dtype=np.bool)
    mask_full = np.zeros((N, N, N), dtype=np.bool)

    full = range(N)
    full.remove(N/2)
    for i in product(full, full, full):
      mask_full[i] = True

    half = range(N/4) + range(3*N/4+1, N)
    for i in product(half, half, half):
      mask_half[i] = True

  else:
    raise ValueError('dimension must be 1, 2, or 3')

  return mask_full, mask_half

Overview¶

PyPFASST is an object-orientated implementation of the PFASST algorithm.

User interactions with PyPFASST are typically marshaled through the PFASST class. This class acts as the overall controller of the algorithm. Implementing a PDE solver in PyPFASST generally consists of the following steps:

Instantiate the PFASST controller.
Add each level of the grid/solver hierarchy, along with appropriate interpolation and restriction operations, to the controller.
Call the run() method of the controller.

Note that if you add only level, then the PFASST algorithm reduces to a serial SDC integrator. In this case it does not make sense to use multiple time processors. If you add multiple levels but only one time processors, then the PFASST algorithm reduces to a space/time multigrid version of SDC (MGSDC).

Levels¶

Each level of the grid hierarchy consists of:

A SDC itegrator. PyPFASST includes once pre-packaged general purpose IMEX SDC integrator:
- IMEXSDC for implicit/explicit schemes.
Users are free to provide their own SDC integrators as extensions of the SDC class.
A function evaluator. These are loosely coupled to the SDC integrator used for the level. For the pre-packaged SDC integrators mentioned above, the user should extend the IMEXFEval to provide their function evaluations.
Spatial interpolation and restriction operators.

Levels are added to the controller from finest (level 0) to coarsest using the controllers add_level() method. Internally, each level is represented by an instance of the Level class.

SDC integrators¶

Each level has an associated SDC integrator. The SDC integrators are implemented as extensions of the base SDC class. User implementations must override the sweep() method.

The constructor of the base SDC class uses the quadrature module to load precomputed SDC integration matrices. The base class also provides a residual() method for computing residuals.

Function evaluators¶

Each level has an associated function evaluator. The function evaluators are implemented as extensions of the base FEval class, but are typically implemented as extensions of either the IMEXFEval class or as dictated by the SDC integrator.

Each function evaluation class must set its shape and size attributes appropriately (see FEval).

Interpolation and restriction operators¶

XXX: description of how the interpolation and restriction routines are called.

Time interpolation and restriction¶

XXX: description of how time interpolation and restriction is done.

Runtime hooks¶

XXX: description of hooks, state etc.

Reference¶

Throughout the PyPFASST documentation, F and G refer to fine and coarse levels.

Top-level classes¶

PFASST class¶

FEval class¶

SDC class¶

Level class¶

Helpers¶

IMEX¶

Low-level routines (guts)¶

Interpolate¶

Restrict¶

State¶

FAS¶

Contributing¶

PyPFASST is released under the simplified BSD license, and is free for all uses as long as the copyright notices remain intact.

If you would like to contribute to the development of PyPFASST, please, dive right in by visiting the PyPFASST project page.