Gallery & Demo Codes

Multigrid solvers: BoxMG in Peano

Together with Marion Weinzierl, I've submitted a paper

 Marion Weinzierl, Tobias Weinzierl: Quasi-matrix-free hybrid multigrid on dynamically adaptive Cartesian grids. ACM Transactions on Mathematical Software

A preprint is available from arXiv. In the paper, we discuss various realisation idioms for multigrid that fuse ideas of FAS, hierarchical generating systems, Galerkin coarse grid operator computations, stencil compression and BoxMG determining problem-dependent inter-grid transfer operators. It should be very straightforward to reconstruct the ideas from the paper at hands of a Peano download. However, I summarise the most important points/pitfalls here:

  • The purely geometric multigrid solvers basically follow the Chapter 4.2 from the Peano cookbook.
  • For the Galerkin and BoxMG part of the manuscript, we add an integer index to each vertex and make it point to a PlainDoubleHeap entry. The stencils are stored on Peano's heap. See the class peano::heap::DoubleHeap for a documentation.
  • To compute Galerkin operators and the BoxMG operators from stencils stored within the vertex can be tedious and error-prone. A similar reasoning holds for the computation of the stencils. Peano offers a toolbox matrixfree that provides routines for this. See the download section or sourceforge. Notably, consult
    • matrixfree::stencil::StencilFactory for routines that yield standard FEM stencils for diffusion and convection operators.
    • matrixfree::stencil::ElementMatrix for routines that allow you to extract from the 4/8 vertices of a cell the element stiffness matrices.
    • matrixfree::solver::Multigrid for routines that accept prolongation, restriction and operator stencils and return the Galerkin coarse grid operators. The class also provides dummies that make prolongation and restriction be d-linear operators.
    • matrixfree::solver::BoxMG offers routines that accept the stencils of one fine grid patch as one large vector and return the BoxMG operators (or simple aggregation or injection).
    The classes and functions in the namespace matrixfree have been built with Peano in mind but they do not directly use Peano for most routines. Instead, they rely on simple double vectors as input. As a result you can also use them in different projects.

Helmholtz problems: geometric additive multigrid, hierarchical basis and BPX-type solvers

This work is based upon work I've submitted together with Bram Reps as

   Bram Reps, Tobias Weinzierl: Complex additive geometric multilevel solvers for Helmholtz 
   equations on spacetrees. ACM Transactions on Mathematical Software

You can read a preprint at arXiv. A source code snapshot of the code used for the paper is available here. Please note that this is a complete snapshot of Peano with all of its glue code and the matrixfree toolbox which is a small collection of routines I use for matrix free linear algebra. As it is a snapshot, both Peano and its toolbox might be available as newer versions today, and it might be reasonable to rewrite a solver from scratch if you need one for a particular project. This snapshot is merely made available as a documentation and to facilitate reproducability. Here are some things to do if you work with the snapshot:

  1. Download the code snapshot from Peano's demo code section.
  2. Untar the archive with
    tar -xzvf helmholtz-snapshot.tar.gz
  3. Type in
    which gives you a list of available build targets. You might have to update the makefile however to adopt pathes or variables. For a quick start with the GNU compiler, try
    make gcc-release-dim2
  4. There is a script create-main-executables that builds (with Intel) all the executables used in the paper. There are release and asserts variants, we have dimensions 2 through 7, there are vectorised (with pragmas) and not vectorised versions, and we have multichannel variants where the number of channels is set through the environment variable CHANNELS. Not all possible combinations are currently provided through the makefile, but the extension of the file should be straightforward if you need additional targets.
  5. Once you have an executable at hand, run the code
    without any argument. It displays a usage message.
    ./peano-2d-gcc-release help
    yields more verbose output.
  6. A first experiment might be
    ./peano-2d-gcc-release sin transition 0.7 0.1 0.01 yes no 100 101
    Our code writes binary legacy VTK files that you can open with Paraview or VisIt, e.g. The file sequence of interest in the first place is channel-0-finegrid*. The files hold one snapshot per iteration, i.e. you can investigate the evolution of the adaptivity pattern.
  7. If you want to study the algorithm's implementation, the most interesting file should be helmholtz/mappings/AdditiveMultigridWithJacobi.cpp and variants of this file.

Compatibility remarks: We use some C++14 features in this code. Our tests with GCC 4.8.3 and GCC 4.9.1 have been successful, but we know that older compiler versions (anthing older than 4.7) might face problems. So please ensure that you have an up-to-date translator version.

Particle in Cell methods

Illustrations stem from the papers

   Weinzierl, T., Verleye, B., Henri, P. & Roose, D. (2016). Two Particle-in-Grid Realisations on Spacetrees. 
   Parallel Computing 52: 42-64.

   Eckhardt, W., Glas, R., Korzh, D., Wallner, S. & Weinzierl, T. (2016), On-the-fly memory compression for 
   multibody algorithms, in Joubert, G.R., Leather, H., Parsons, M., Peters, F. & Sawyer, M. eds, Advances 
   in Parallel Computing 27: International Conference on Parallel Computing (ParCo) 2015. Edinburgh, Scotland, 
   IOS Press, Amsterdam, 421-430.

You can find a prepring of the Parallel Computing manuscript at arXiv. The underlying particle management within Peano is available as predefined user toolbox: if you want to use it, its merely a few simple additions to the specification file. No particle update has to be rewritten.

PeanoClaw/hyperbolic solvers

Illustrations stem from the dissertation

 K. Unterweger: High-Performance Coupling of Dynamically Adaptive Grids and Hyperbolic Equation Systems, 2016
The underlying patch management within Peano is available as predefined user toolbox: if you want to use it, its merely a few simple additions to the specification file. This toolbox is also used in the ExaHyPE project. The shockbubble experiment above uses the 3D Euler solver of David Ketcheson as found in the Clawpack solver. The tsunami simulation uses shallow water kernels by Michael Bader's group with a bathymetry from GEBCO ( The GEBCO_08 Grid, version 20100927, ).

Shallow water/experimental code

The video discusses work of the paper

 Weinzierl, Tobias, Bader, Michael, Unterweger, Kristof & Wittmann, Roland (2014). Block Fusion on Dynamically Adaptive Spacetree Grids for Shallow Water Waves. Parallel Processing Letters 24(3): 1441006.