SummationByPartsOperators.jl: A Julia library of provably stable discretization techniques with mimetic properties

SummationByPartsOperators.jl is a Julia library of summation-by-parts (SBP) operators, which are discrete derivative operators developed to get provably stable (semi-) discretizations, paying special attention to boundary conditions. Discretizations included in this framework are finite difference, Fourier pseudospectral, continuous Galerkin, and discontinuous Galerkin methods. The main aim of SummationByPartsOperators.jl is to be useful for both students learning the basic concepts and researchers developing new numerical algorithms based on SBP operators. Therefore, SummationByPartsOperators.jl provides a unified framework of all of these seemingly different discretizations. At the same time, the implementation is reasonably optimized to achieve good performance without sacrificing flexibility.


Statement of need
Partial differential equations (PDEs) are widely used in science and engineering to create mathematical models of real-world processes. Since PDEs often need to be solved numerically, a vast amount of numerical methods has been developed. Since it is impossible to keep up with all recent research, sub-communities focussing on specific applications or numerical methods emerged. Sometimes, these communities develop different vocabulary and notations, making it hard for newcomers (or even experienced researchers) to see similarities and connections between seemingly unrelated approaches. To transfer new ideas and developments and knowledge from one community to another, common abstractions can be helpful. The concept of SBP operators is such an abstraction. In recent years, SBP operators have attracted a lot of attention, in particular for PDEs modeling advection-dominated problems, where they enabled the construction of energy-or entropy-stable numerical methods, including finite difference, discontinuous Galerkin, continuous Galerkin, and (pseudo-) spectral methods. Their success is based on mimetic properties which enable the transfer of results obtained for differential equations at the continuous level to the discrete level. In particular, SBP operators are designed to mimic integration-by-parts discretely as summation-by-parts, enabling discrete analogues of energy/entropy methods for PDEs.
SummationByPartsOperators.jl is written entirely in Julia (Bezanson et al., 2017). Making use of multiple dispatch and generic types, SummationByPartsOperators.jl provides a unified interface for different SBP operators. At the same time, the implementations are reasonably fast (again, due to multiple dispatch and specialized implementations for each operator class). Together, this facilitates the development of new algorithms and research in numerical analysis, which is the primary goal of this package. In addition, SummationByPartsOperators.jl has been used in a number of graduate-level numerical analysis courses, allowing students to understand the connections between different SBP methods by presenting them in a unified framework. In addition, some of the operators were not available in open source software previously (to the best of the author's knowledge).

Features
SummationByPartsOperators.jl implements numerical methods based on SBP operators of the following classes: • finite difference methods • Fourier collocation methods • continuous Galerkin methods • discontinuous Galerkin methods Since a discrete derivative operator is a linear operator, all of these SBP operators implement the basic interface of such linear operators (AbstractMatrix in Julia) such as multiplication by vectors and addition of operators. Finite difference and Fourier operators on periodic domains also allow the construction of rational functions of operators and their efficient implementation using the fast Fourier transform (Frigo & Johnson, 2005).
In addition to basic SBP derivative operators, SummationByPartsOperators.jl contains a number of related operators, such as • SBP artificial dissipation operators • spectral viscosity operators for Fourier methods • modal filter operators for Fourier methods and Legendre pseudospectral methods Using Julia's rich type system, all of these operators are implemented as their own types. This enables several optimizations such as a memory requirement independent of the number of grid points. In contrast, implementations based on sparse/banded matrices have a memory requirement growing linearly with the number of grid points. In addition, the structure of the operators can be taken into account for operator-vector multiplications, usually resulting in speed-ups of an order of magnitude or more on consumer hardware. For example, the application an optimized the sixth-order (in the interior) finite difference SBP operator (Almquist, 2017) on a grid with 1000 nodes takes roughly 330 ns on a consumer CPU from 2017 (Intel® Core™ i7-8700K) using version v0.5.5 of SummationByPartsOperators.jl. In contrast, the same operation takes roughly 3.9 microseconds using a sparse matrix format used in other implementations of this operator (Almquist, 2017). This benchmark is based on the following code, which also provides a very basic example of SummationByPartsOperators.jl. The output of the last command will be a relatively small number on the order of 4.2e-13.
Following good software development practices, SummationByPartsOperators.jl makes use of continuous integration and automated tests required before merging pull requests. Documentation is provided in form of docstrings, a general introduction, and tutorials. In addition, SummationByPartsOperators.jl is a registered Julia package and can be installed using the built-in package manager, handling dependencies and version requirements automatically.

Related research and software
There are of course many open-source software packages providing discretizations of differential equations. However, many of them focus on a single class of numerical methods or a specific application, e.g., • finite difference methods (DiffEqOperators.jl, a part of DifferentialEquations.jl (Rackauckas & Nie, 2017)) • finite volume methods (Oceananigans.jl (Ramadhan et al., 2020), Kinetic.jl (Xiao, 2021)) • spectral methods (ApproxFun.jl (Olver & Townsend, 2014), FourierFlows.jl (Constantinou & Wagner, 2021)) • finite element methods (Gridap.jl (Badia & Verdugo, 2020)) • discontinuous spectral element methods (Trixi.jl Schlottke-Lakemper et al., 2021 We are not aware of any open-source software library implementing all of the SBP classes using a unified interface or even several finite difference SBP operators on finite domains, which are usually heavily optimized (Mattsson et al., 2014(Mattsson et al., , 2018 and not available in other open source packages. Sometimes, restricted sets of coefficients are available online (Almquist, 2017;O'Reilly, 2019), but there is no other extensive collection of these methods.
Of course, there is a plethora of additional open source software implementing numerical methods for PDEs and each package has its own design criteria and goals. SummationBy-PartsOperators.jl provides a unified interface of different SBP operators. Thus, there is a partial overlap with some of the packages mentioned above such as finite difference operators on periodic domains (DiffEqOperators.jl, but with a different handling of bounded domains) or Fourier methods on periodic domains (ApproxFun.jl, but with a different interface and extensions). In addition, many packages focus on a specific application such as some specific fluid models (Oceananigans.jl, FourierFlows.jl) or hyperbolic PDEs (Trixi.jl). In contrast, SummationByPartsOperators.jl focuses on the numerical methods and provides them in a form usable for rather general PDEs. For example, there is ongoing work to use the basic operators provided by SummationByPartsOperators.jl in Trixi.jl.