pystencils+lbmpy: An Algebraic Language for Stencil-Based Numerical Methods

Stencil schemes are a family of numerical methods where the values of mathematical functions (fields) on a computational grid are updated via equations involving only values from neighboring grid primitives. Among them count classical finite-(difference/volume/element) approaches to solving PDEs, as well as the lattice Boltzmann method (LBM). The efficient evaluation of stencil operators is thus a critical building block for many numerical solvers. This lecture introduces the software packages pystencils and lbmpy, which together offer an algebraic toolbox for the design of stencil schemes in general, and lattice Boltzmann methods in particular. It will introduce the fundamental concepts of both packages, which attendees will later interactively explore, and apply to actual numerical problems, during the hands-on sessions. The centerpiece of pystencils is a purely algebraic language where grids, fields, and a stencil's equations are represented as symbolic objects. Based on the computer algebra system SymPy, users write numerical operations on a mathematical level of abstraction. Pystencils assists this process by offering auto-discretization for finite-difference and finite-volume schemes. Lbmpy extends the language with concepts specific to the LBM, and exposes interfaces for the definition of lattice Boltzmann methods by their mathematical components. Central to lbmpy is its algebraic engine, by which implementations of the key LBM operations (initialization, collide-stream, forcing, boundary conditions, etc) are automatically derived. It supports a wide range of LBM schemes, augmented by various boundary handling and forcing models. The algebraic operators written by the user and/or derived by lbmpy are finally passed to pystencils' code generation system, which transforms them into executable code targeting modern parallel CPUs and GPU accelerators. The code generator exploits available domain-specific knowledge to maximize the kernels' runtime efficiency on the available hardware, using transformations such as memory-friendly loop ordering; use of OpenMP and CPU vector intrinsics; and hierarchical reductions on GPU. Generated operator code can finally be employed to run simulations entirely within Python, where Pystencils supports the rapid prototyping of numerical methods through a runtime system based on NumPy and CuPy. Kernels may also be exported for use with external C++ code bases (such as the waLBerla framework) in order to build larger, production-grade simulation applications.