 |
Peano
|
|
Peano
|
All numbers stem from an eight-core workstation, the profile variant (domain 54, h=1, t=5), eight OpenMP threads, and a Peano build as described above. "Sweeps" is the number of mesh traversals of the time stepping, i.e. the number of time steps (the time step size follows the CFL condition).
| Scheme | Wall time [s] | Sweeps | Time step size dt |
|---|---|---|---|
| fv-regular | 9.3 | 35 | 0.15 |
| fv-amr-1 (radius 1) | 174 | 102 | 0.05 |
| fd4-rk2 | 31 | 64 | 0.083 |
| fd4-rk3 | 47 | 64 | 0.083 |
| fd4-rk4 | 62 | 64 | 0.083 |
| fd4-rk4-lopsided | 82 | 64 | 0.083 |
| fd4-rk4-lopsided-ko | 82 | 62 | 0.083 |
(fv-amr-1 and the two lop-sided runs use --compiler sequential, see the warning below. The other runs use the default compiler producer-consumer-fusion.)
What we learn from this:
--variant production --end-time 100 (and checkpoints) on a cluster for that.producer-consumer-fusion if more than one thread is used (all threads of the runs with 2, 4 and 8 threads stood still without a single time step; with one thread the same executable makes progress). The sequential compiler (--compiler sequential) runs fine. See the profiling section below for the diagnosis. The first run of the scheme study consequently hit the timeout for these two schemes. The runs do not write checkpoints. To compare solution profiles in addition (lapse over radius), pass --extra-arguments "--checkpoint-interval 1" and extract the profiles with probe-single-black-hole.py.
Fourth-order Finite Differences with AMR will be added to this study once the corresponding branch is merged.
runs the Finite Volumes on the regular mesh and with one AMR level (regular, amr-1; radius 1) for each task graph compiler (sequential, fusion-only, producer-consumer-sequential, producer-consumer-fusion) and every thread count. The thread count is set via OMP_NUM_THREADS and therefore applies to OpenMP builds only. The profile variant uses Peano's stats build. The study measures the wall time of each run (results-performance.csv).
Things to compare:
Wall time of the whole run in seconds (setup, i.e. TwoPunctures and mesh construction, plus 35 time steps on the regular mesh or 102 on the adaptive mesh), eight-core workstation, Finite Volumes, profile variant. The regular mesh figures are the mean of two runs (they differ by less than 0.5 s), the adaptive mesh figures stem from one run.
| Mesh | Compiler | 1 thread | 4 threads | 8 threads | Speedup 1 to 8 |
|---|---|---|---|---|---|
| regular | sequential | 39.0 | 14.5 | 12.1 | 3.2 |
| regular | fusion-only | 43.8 | 13.4 | 9.4 | 4.7 |
| regular | producer-consumer-sequential | 39.4 | 14.5 | 12.4 | 3.2 |
| regular | producer-consumer-fusion | 45.3 | 12.2 | 9.5 | 4.8 |
| adaptive (1 level, radius 1) | sequential | 511 | 172 | 3.0 | |
| adaptive (1 level, radius 1) | fusion-only | 584 | 130 | 4.5 | |
| adaptive (1 level, radius 1) | producer-consumer-sequential | 504 | 147 | 3.4 | |
| adaptive (1 level, radius 1) | producer-consumer-fusion | 592 | 105 | 5.6 |
producer-consumer-fusion) is the fastest as soon as more than a handful of threads are used: 1.3 times faster than the sequential compiler on the regular mesh and 1.6 times faster on the adaptive mesh with eight threads. It is, however, the slowest variant with one thread (15% slower on both meshes), i.e. it pays an overhead that only parallelism recovers.producer-consumer-sequential) is no faster than the sequential compiler on the regular mesh (12.4 s vs. 12.1 s) and only 17% faster on the adaptive mesh. The gain comes from fusing sweeps, not from tasking: fusion-only (no tasks) matches producer-consumer-fusion on the regular mesh (9.4 s vs. 9.5 s). Tasks gain another 24% on the adaptive mesh (105 s vs. 130 s), where the patches are costlier and the tasks have more to overlap.producer-consumer-fusion.For the compute-bound view of the kernels, see Benchmarks of Compute Kernels of CCZ4.
We profiled the profile variant with Intel VTune (2024.3, hotspot analysis with software sampling, i.e. no hardware counters are needed) on the eight-core workstation. The Peano build is the one described above, with -g -O3. A profile looks like
The percentages below refer to CPU time summed over all threads. We profiled five setups:
| Setup | CPU time [s] | Remarks |
|---|---|---|
FV, regular mesh, producer-consumer-fusion, 1 thread | 47 | cleanest view on the kernels |
FV, regular mesh, producer-consumer-fusion, 8 threads | 66 | spin time 3.1 s, overhead 3.5 s |
FV, regular mesh, sequential, 8 threads | 86 | spin time 6.9 s, overhead 7.5 s |
FV, adaptive mesh (1 level), producer-consumer-fusion, 8 threads | 837 | spin 40 s, overhead 47 s |
FD4 (RK4), regular mesh, producer-consumer-fusion, 8 threads | 477 | spin 35 s, overhead 38 s |
Finite Volumes, one thread.** The three CCZ4 routines plus the generated Rusanov kernel make up 63% of the time: exagrype::ccz4::ncp (24%), the generated kernel CCZ4_timeStepWithRusanov itself (23%, the ncp implementation is partially inlined into it), exagrype::ccz4::source (12%), nonconservativeProduct (3%). TwoPunctures (PunctEvalAtArbitPositionFaster, the initial data) takes 18% of this short run. Heap management takes 8%: free 4.3%, operator new 2.2%, memset 1.2%. memcpy takes 2.4%. Peano's mesh traversal (peano4::grid) is below 2%.
Finite Volumes, eight threads.** The compute kernels shrink to 39% (they scale), whereas everything that moves data does not: memcpy grows from 1.2 s to 5.1 s of CPU time (8%) and peano4::grid takes 3.1 s with producer-consumer-fusion but 7.9 s (9%) with the sequential compiler. That is the effect that makes fusion fast: it cuts the traversal cost by 60%. The remaining 6 s of spinning at barriers (__kmp_fork_barrier) and the serial setup limit the speedup.
Adaptive mesh.** Compute kernels dominate even more (CCZ4_timeStepWithRusanov 17%, source 13%, nonconservativeProduct 4%), but memcpy (8%), TwoPunctures (5%, the new volumes need their initial data), getVertexType (4%) and allocation (free, memset, new, memalign: together 3%) still take a combined 20%. The setup takes 8.7 s of the 105 s.
Fourth-order Finite Differences.** memcpy is the biggest single item with 14% (65 s), followed by ncp (13%), source (10%), peano4::grid::GridTraversalEventGenerator::getVertexType (6%, an event generator which should cost nothing) and the right-hand side kernel itself (4%). Spin and overhead time account for another 15% (73 s).
The Finite Differences with lop-sided advection (--lopsided) deadlock with the compilers that spawn tasks as soon as OpenMP uses more than one thread. We attached gdb to a hanging run with eight threads: all eight threads sit inside tarch::multicore::omp::internal::hasTaskFinishedInOpenMP5() (src/tarch/multicore/omp/Tasks.cpp), called from the body of a task created in spawnTaskWithOpenMP5(). That function polls the task's input dependencies in a busy loop inside the task body without ever yielding. If every thread executes a task that waits for a predecessor which still sits in the task queue, nobody is left to run the predecessor. The lop-sided advection adds an action set AddLopsidedAdvection_WaitForTask to the task graph (the plain FD4 scheme has none), which presumably creates the dependent chains that trigger this situation; we did not analyse the dependency structure further. With one thread, or with the sequential compiler, the same code runs through. We consider this a bug of the OpenMP 5 fallback (the compiler offers OpenMP 5, so the OpenMP 6 dependency path is not used).
Sorted by expected benefit per effort; we did not change any code in Peano.
#pragma omp taskyield to the loop that polls hasTaskFinishedInOpenMP5() in spawnTaskWithOpenMP5() (the polling loop for the in/out dependency a few lines above it already has one) or defer the task until its predecessors have finished. Until then, run the Finite Differences with lop-sided advection with --compiler sequential. We did not test the taskyield patch.producer-consumer-fusion for runs with at least four threads and sequential for a single thread (see the table above). The default of the script is producer-consumer-fusion, which is the right choice for the production runs but not for debugging on a single core.kernels/CCZ4_rusanov3d.cpph) allocates 30 temporary arrays per call, each as an array of pointers plus one new double[...] per cell of the batch, and deletes them again. For CCZ4 with 59 unknowns, this costs 8% of the single-thread time (plus the memsets that initialise them). Reusing thread-local scratch buffers, or fusing the loops so that the lambda and flux fields stay in registers, removes that. It lives in RusanovWithExplicitEuler.py and the DSL back end, i.e. in the code generator and not in the benchmark.ncp (24%) and source (12%) are evaluated three times (once per face direction) on averaged states, and maxEigenvalue is evaluated six times per volume (three times on the halo extent for the Rusanov flux and three times on the result for the next time step). The latter is a cheap routine for CCZ4, but the former two are the most expensive in the profile. Check if their arithmetic vectorises (use the optimisation report as described in Benchmark implementation); we did not run a micro-architecture analysis.memcpy takes 8% on eight threads for the Finite Volumes and 14% for the Finite Differences, and it is the one item that scales worst with the number of threads (it is probably limited by memory bandwidth). The call stacks show that both for the Finite Volumes and for the Finite Differences, the copies come from toolbox::blockstructured::api::ReconstructPatchWithHalo (its touchCellFirstTime()): in every time step (every Runge-Kutta stage for FD4), each cell copies its patch and the halo of its neighbours into a temporary array before the kernel runs. Reading the halo straight from the face data inside the kernel, or reusing the reconstructed patch between the stages, would save bandwidth. We did not look into how the copy is implemented in detail.getVertexType in the Finite Difference profile (6% of 477 CPU seconds for a function that classifies vertices): with eight threads it competes with the kernels; presumably it is called once per vertex and per Runge-Kutta stage. A cached classification would remove it. We did not verify the call count.The VTune results are not part of the repository. The command lines above reproduce them.
One difficulty in letting ExaGRyPE run is that there's a plethora of different setting that one can pick from. This is quite natural for a research code but quickly becomes challenging for people who just wanna have a go or start benchmarking. It also means that shipping a proper benchmark is far from trivial.
The single black hole benchmark is written in a way such that it supports many different settings. We can therefore play around with them and find good choices. Anecdotal evidence fees into the settings used by benchmarks_exahype2_ccz4_single_black_hole_production_benchmarks. All the different benchmarking ingredients for the HPC assessment are available through the script performance-studies.py.
This page first provides some generic description of the variants. This text serves primarily as background reading. Most readers might skip it. From hereon, we demonstrate how to assess certain system properties systematically:
Each section discussing one of the items above kicks off with a reference configuration. You should be able to run this as it is as kind of a sanity check.
Benchmarking and playing around with the single black hole is never easy, as it is a rather complex setup already with a rather long initialisation and grid construction phase. However, it is possible to benchmark it once you disable plots, the tracers, and you use an unreasonable coarse mesh. That is, we disable a lot of key features for the benchmarking runs. In general, we apply the following simplifications:
-et 1e-2 via the Python API front-end.--plot 0 from the Python call to zero which switches off the plotting.A typical setup therefore resembles
The script allows us to toggle between different solver variants from the command line by picking a variant through --solver.
fv- prefix (if you use the command line) followed by the patch size. We exclusively employ enclave solvers in this setup, i.e. each Finite Volume solver spawns tasks.Proper load balancing between ranks and threads is obviously key for good performance. The benchmark supports a couple of load balancing schemes to empirically study different concepts:
generic uses an instance of toolbox::loadbalancing::strategies::cascade::SpreadOutHierarchically. That is, we first populate all ranks and then all threads, unless the Python API specifies the value trees which constraints the maximum number of subpartitions (trees) that we employ per rank.tailored is generic as well, but it uses a cost functional which hard-codes information about the Schwarzschild solution: The closer a cell to the centre, the higher its weight.offline-regular uses the offline balancer over a regular mesh, i.e. it uses the cell size to anticipate what a regular mesh would look like and then computes an optimal partitioning before we even start the compute.Constraining the number of trees does not mean that threads idle: We use task-based parallelism on top of domain decomposition and hence assume that idle threads grab tasks and hence automatically balance out illpartitioning.
For each benchmarks, there are various schedulers:
--scheduler native Picks the native, hardcoded scheduler from Peano. This disables all task fusion and also all GPU offloading.--scheduler native-piorities Constructs an instance of benchmarks::exahype2::ccz4::MulticoreOrchestration as task orchestration. It still disables all the task fusion (and therefore GPU offloading), but now activates features such as--scheduler gpu-lazy-evaluation Tasks are shipped off to the GPU lazily, i.e.~when the task results are requested. Compare to the description in tarch::multicore::orchestration::FusionImplementation.--scheduler gpu-timed-execution Tasks are shipped en bloc to the GPU after each mesh sweep that produces tasks. Compare to the description in tarch::multicore::orchestration::FusionImplementation.--scheduler gpu-producer-consumer Shipped off to the GPU with producer-consumer tasksCompare to the description in tarch::multicore::orchestration::FusionImplementation.--scheduler gpu-eager-consumer-tasks Tasks are shipped off to the GPU using Peano's eager consumer pattern.Compare to the description in tarch::multicore::orchestration::FusionImplementation.--scheduler --no-tasks No tasking is used, i.e. everything that could be mapped onto a task is executed immediately.If an orchestration with task fusion is picked and you employ a combined FV-FD solver, then only the Finite Volume tasks are sent off to the accelerator. Details on this behaviour can be found in benchmarks::exahype2::ccz4::MulticoreOrchestration.
We start with some experiments where we artifically restrict the number subpartitions uses on the node. If we use one subpartition per core, we quickly run into a situation where the individual subpartitions become too small and overhead eats up all scalability improvements. Therefore, a first step restricts the number of trees through --trees. While we do so, we first disable all tasking manually, and we also stick exclusively to Finite Volumes:
There are two job script in the scripts subdirectory, which can serve as a starting point if you want to submit experiments to a supercomputer with an array job. This script is for Hamilton 8 but most Slurm scripts will look similar.
Once the data is in, it might be best to package it into two archives
and to copy them over onto a local computer. There is a dedicated script plot-trees-vs-runtime.py in the scripts subdirectory of the benchmark that can read zips and create output data. The script is based upon the load balancing postprocessing scripts and yields plots similar to the one below.
The above images stem from an OpenMP run on Hamilton 8 (top) and a run on the same machine using Intel's TBB with Peano's TBB extension (bottom).
We see that there is a sweet spot whenever we hit the canonical \( 3^d \) subpartition count. We also see that there is limited sense in going beyond 27 subpartitions on a rank. The crucial question here is whether this is a general pattern, or if the workload per partition then simply becomes too small.
To a certain degree, we have to analyse adaptive meshes to find an answer to this question, although an in-depth performance analysis would obviously be even better.
The above images stem from an OpenMP run on Hamilton 8 (top) and a run on the same machine using Intel's TBB with Peano's TBB extension (bottom).
The adaptive mesh data is more interesting: The balancing of cells with the present load balancing scheme is poor. Consequently, we need a more sophisticated load balancing scheme.
Once we are reasonably satisfied with the load balancing per rank, it is an interesting question what kind of tasking back-end is most advantageous for the settings. That is, now we fix the number of subpartitions to a reasonable choice and we alter the scheduling orchestration. We know that we cannot exploit all cores through domain decomposition alone, and hence hope that the tasks step in and help us exploit those cores that we haven't used yet, i.e. we use the whole node, use a reasonable number of subtrees per node and hence use a few cores to handle those subtrees, and then let those subtrees spawn tasks which should bring in the idling remaining threads.
Despite the prefix gpu, all of these schemes also should work on a CPU-only system. Again, there's an example Slurm submission script in the scripts subdirectory which works for Hamilton. This time, we cannot straightforwardly use array jobs, so I decided to stick with one big script which runs through the variants. We do not farm out.
The mesh size is delibertely chosen coarser than in a previous setup. As we bundle tasks, we temporarily need more memory, so having some slack would help.
Again, the cleanest way is to create on archive as my postprocessing script can handle archives directly:
This time, I pack the regular grid data and the adaptive grid data into one single archive. The script should visualise both variants, and it also provides a plot with a linear y-axis and one where we use the log scale.
| Scale | OpenMP | TBB |
|---|---|---|
| linear | 
| 
|
| log | 
| 
|
The above images stem from an OpenMP run on Hamilton 8 (first two plots) and a run on the same machine using Intel's TBB with Peano's TBB extension (bottom two plots).
The lessons learned now are not trivial anymore:
Below is a template of a script how to creates data points. You will have to amend the prefix accordingly. Note that we really only are interested in one line from the output and that's the time to solution that is spent on the actual solver, i.e. bare any I/O and mesh generation.
The last line packs the data from one upscaling exercise into a tar archive. Single files are obviously more convenient to move around between computers, but archives are notably directly supported by Peano's postprocessing scripts.
Get the tar archives. After that, we can use the out-of-the-box Peano scripts to generate the plots
which clearly showcase that the kernel parallelisation pays off even though the overall scalability is poor with the present rather coarse mesh:
Here are some canonical settings:
| Description | Solver | Cell size | End time | Patch size | Remarks | Build command |
|---|---|---|---|---|---|---|
| Minimal setup | Finite Volumes | 5.4 | 1.0 | 6 | Toy problem with 729 cells. Tells you if the code, in principle, works. Should crash at one point suddenly as the eigenvalue becomes zero, which showcases that this setup is too coarse to capture anything that has a physical meaning. | python3 performance-studies.py -s fv-6 -et 1.0 --plot 0.0 --cell-size 5.4 |
| Simplest FV setup for single node workstation (laptop) | Finite Volumes | 1.8 | 1e-1 | 6 | You cannot see AMR here, but already gives you 19,683 cells. Load balancing should be really good. This is kind of the largest regular mesh setup that I can fit on a local laptop. | python3 performance-studies.py -s fv-6 -et 0.1 --plot 0.0 --cell-size 1.8 |
1.9.8