Sparse solver benchmarks

For selection and tuning of specific sparse solvers, it is good to have relevant benchmarks. As indicated in ‘Solving PDEs’, the performance of a solver depends on the type of equation being solved, the mesh geometry, boundary conditions. Different solver algorithms have different strengths and weaknesses, and should be tested for different transport phenomena.

Here, we collect relevant benchmarks, i.e. actually solving PDEs appearing in the study of transport phenomena.

PyFVTool sparse solver benchmark 260826

Comparison between built-in solver and external solver. Simple wall time measurement.

• SuperLU (single-threaded): scipy.sparse default built-in spsolve (SciPy 1.18.0)

• PARDISO (multi-threaded): Intel MKL PARDISO via pypardiso (MKL 2026.0.1)

Hardware: MOLTECH-Anjou HPC compute node (32+ Gb RAM, Intel Xeon, 20+ threads, Linux)

SLURM/environment config: Ncpu=32; RAM=32Gb; MKL=unrestricted

Calculation: time-dependent advection-diffusion on 2D axisymmetric cylindrical grid \((r, z)\)

Grid size	FVM cells	SuperLU	PARDISO	Speed-up
Nr x Nz	incl. ghost	iter/s	iter/s
50 x 2070	107744	2.365	16.90	7.14
50 x 4500	234104	1.118	8.524	7.62
100 x 2070	211344	0.599	8.899	14.86
100 x 4500	459204	0.290	4.060	14.00

SciPy’s built-in scipy.sparse.linalg.spsolve with SuperLU provides a solid, default baseline. Note that it is single-threaded and does not perform any optimization (e.g. re-using symbolic/numerical factorization). Make sure that SciPy actually uses the built-in SuperLU by calling scipy.sparse.linalg.use_solver(useUmfpack=False). Without this call, SciPy’s sparse solving behaviour becomes unpredictable.

Intel’s MKL PARDISO solver was used without explicit control of the number of threads used (Using top, we estimate between 10 and 20 threads, depending on problem size). The parallelization by PARDISO is quite efficient in our 2D cylindrical case.

PyFVTool sparse solver benchmark 260830

Here we keep track of different computations with the MKL PARDISO solver while gradually changing the amount of compute and memory requested. We now explicitly set the number of threads for MKL.

These are simulations for an ongoing research project, and we monitor solver performance as we go via simple wall-time measurements. We also monitor memory, but memory is not really an issue (memory use is modest), so it is not reported here. Furthermore, the calculations are not entirely identical: some simulation parameters vary, but the general FVM computation is very similar. For several runs on the same grid size, we take the smallest value for the compute rate (‘iter/s’), which varies within 5…10%.

Hardware: MOLTECH-Anjou HPC compute node (Intel Xeon)

Calculation: time-dependent advection-diffusion on 2D axisymmetric cylindrical grid \((r, z)\)

Grid size	FVM cells	PARDISO	SLURM/environment
Nr x Nz	incl. ghost	iter/s
75 x 2070	159544	13.01	Ncpu=16; RAM=16Gb; MKL=16threads, dynamic
75 x 4500	346654	5.703	Ncpu=16; RAM=16Gb; MKL=16threads, dynamic
75 x 4500	346654	5.457	Ncpu= 8; RAM= 8Gb; MKL= 8threads, dynamic
75 x 5400	415954	3.787	Ncpu= 8; RAM= 8Gb; MKL= 8threads, dynamic
75 x 4800	369754	6.188	Ncpu=16; RAM= 8Gb; MKL=16threads, dynamic

It seems that restricting MKL PARDISO to 8 threads instead of 16 only has a minor effect on the compute rate for this computation, at similar grid dimensions. This makes it easier to squeeze in the compute job on the HPC cluster. However, 16 threads seems the more generally ‘safe’ option, when the grids are larger. Memory seems not to be an issue: the computations work fine with “only” 8 Gb.

Of note, there is likely some heterogeneity in performance between different compute nodes on the HPC cluster, which we do not control or monitor here. This set of benchmarks is mainly to find reasonable SLURM requests (Ncpu, RAM), while also monitoring operational performance.