Parallelisation and Performance Portability in Schnek
Posted 5th August 2026 by Holger
Some time ago, I reported on my first steps using Kokkos to write code that runs on both GPU and CPU. Kokkos is a library aimed at scientific computing defining multi-dimensional arrays, called views. The library makes it possible to use accelerated hardware, like GPUs, when it’s available, but also allowing to automatically fall back to classic serial execution. The experience I gathered with that library made it an ideal candidate for integration with my Schnek project. The Schnek library also has multidimensional grid and field classes and these classes allow different back-end storage models.
Since then, I have been developing Schnek along two main lines. The first line involved the implementation of Kokkos views as a back-end storage model for Schnek. This also involved supporting a kernel based iteration over a multi-dimensional region of a grid. In future, simulation codes will not contains the iteration loops themselves. They will specify ranges and hand over a kernel function to the iteration primitive which then calls the kernel function for each index in the range. Depending on the execution back-end, this can result in classical sequential iteration or a parallel dispatch on some accelerated hardware.
The second line of development of Schnek involves a re-thinking of the current MPI domain subdivision. The class structure of the subdivision made it difficult to move to more complex MPI parallelisation schemes. This meant implementing a completely new domain decomposition paradigm. Similar to the callback based philosophy of the new iterations, more control is handed over to the new domain decomposition strategy. Simulation code will not be expected to allocate grids of fields itself. Instead, it will supply factory functions to the decomposition strategy and obtain registration handles. It can then create iteration contexts with multiple handles and provide callback functions that will be called for each n-dimensional region.
As part of the new domain decomposition strategies, I already implemented a function to convert indices to positions on a Hilbert curve. This was more of a mental exercise at the time but also a way of manifesting my intention of implementing a parallelisation strategy based on space filling curves. This function is, as of now, unused but I am fully intending to implement the decomposition that uses it.
Both of these tasks are almost done, but yet not fully completed. But one big task is currently still remaining. Joining the two lines of development together. Making the new domain decomposition work with Kokkos is not a trivial undertaking. How do you handle boundary exchange for grids that are stored on GPU memory. Modern MPI versions support device memory exchange. But the library needs to create the exchange buffers on the correct memory space and copy them between the grid and the buffers. Alternatively, for CPU exchange I am experimenting with using strides to avoid copying between the grids and separate memory buffers altogether.
I also started implementing an ambitions declarative algorithm specification, but that will probably not make it into version 2.0 of Schnek.
In the following weeks and months, I intend to post more regularly on the progress of the development of Schnek. My intention is to communicate design decisions and trade-offs as well as ideas for future development.
Leave a Reply