Code structure and parallelization

The package struphy is organized in layers. A Simulation combines a model, a domain, a grid and the numerical options. It builds the discrete de Rham complex on the domain and grid (feec), instantiates the model's variables (finite element coefficient vectors or marker arrays), and then calls the model's propagators in sequence at every time step. A propagator either solves a linear or nonlinear system on the finite element spaces, or calls a pusher or accumulator that applies a compiled kernel to all markers.

Struphy relies on a few external packages for the infrastructure below this layer.feectools, a fork ofPsydac, provides the spline spaces, the distributed stencil vectors and matrices, and the MPI wrapper. Pyccel compiles the kernels. cunumpy provides the array module used throughout the code (NumPy or CuPy, see below), and scope-profiler records timings for the profiling regions.

The kernels are ordinary Python functions with type annotations, restricted to a NumPy subset that Pyccel can translate. There are roughly three dozen kernel modules, in bsplines,feec, geometry, linear_algebra and pic: spline evaluation, domain mappings, mass-matrix assembly, matrix–vector products, marker pushing, and the accumulation of marker contributions onto the grid. Working with a single Python source means that a kernel can be read, tested and, if needed, debugged as Python.

struphy compile (or the Compiler class in Python) transpiles the kernels to Fortran or C and builds shared libraries. Kernel modules that depend on each other are rebuilt in the right order and only when needed. The GNU compiler is the one regularly tested; Intel, PGI, NVIDIA and LLVM compilers can be selected through the --compiler option or a compiler configuration file. A kernel that has not been compiled still runs as pure Python, which is slow but useful for development.

Three mechanisms are available and can be combined. They address different limits: OpenMP uses the cores of one node, domain decomposition distributes both the fields and the markers, and clones distribute only the markers when the field grid does not need to be split further.

The logical grid is divided into blocks, one per rank, and the finite element vectors and matrices are stored in distributed stencil format with ghost regions at the block boundaries. The directions that are decomposed are chosen with mpi_dims_mask. As a guideline, a direction should only be decomposed if every rank keeps enough elements in it.

Each marker belongs to the rank whose block contains it. After the markers have been moved, they are sorted according to the decomposition (mpi_sort_markers): markers that left a block are sent to their new owner, taking periodic boundary conditions into account. Accumulation kernels deposit marker contributions into the local block and its ghost region, followed by a ghost-region exchange. For SPH models, ghost markers are exchanged between neighbouring blocks.

With num_clones greater than one, the MPI ranks are split into groups. Each clone carries the field decomposition and its share of the markers, and the accumulated quantities are summed over the clones (an Allreduce over the communicator that links equivalent ranks of different clones). The number of ranks must be divisible by the number of clones. This increases the number of markers that can be used without splitting the field grid further, at the cost of replicating the fields.

Markers are stored in a single two-dimensional array, one row per marker, with columns for position, velocity, weights and buffers used by the pushers. Removed markers leave holes that are refilled when markers arrive from other ranks. With the option sort_step, markers are also sorted in memory every N steps so that markers in the same cell are close together, which improves cache use in evaluation and accumulation.

A run can record the time spent in the setup, in each propagator, pusher and accumulation kernel, in marker sorting and communication, and in the linear solves. The results are written to an HDF5 file, from which per-rank timelines (Gantt charts) and flame graphs are generated. Optional hardware-counter (LIKWID) and NVIDIA (NVTX, CUDA events) instrumentation is available in the profiler options. The examples on this site show these charts for the example runs.

Building a simulation lists the profiling regions and the options of a run.

Current status. The array module xp from cunumpy selects NumPy or CuPy as the array backend, so that the Python part of the code can allocate its arrays on a GPU. The compiled kernels, however, still run on the CPU. Pyccel produces Fortran or C code that only accepts NumPy arrays, so CuPy arrays are copied to the host before a kernel call and the results are copied back afterwards. The CuPy backend therefore does not accelerate the kernels by itself, and the numerical results do not depend on the backend. The profiler can already record NVTX ranges and CUDA-event timings.

Planned. The particle kernels (pushers and accumulation) dominate the cost of kinetic simulations, and porting them to run on the device is the next step. The approach under consideration is to implement them as CuPy RawKernels (CUDA C), to be used in place of the Pyccel kernels when the CuPy backend is active, so that markers and fields stay on the device between substeps.