Pyroclast Benchmark Eval

---
name: pyroclast-benchmark-eval
description: Evaluates the numerical correctness, stability, and parallel performance (CPU/GPU strong scaling and distributed weak scaling) of a matrix-free finite difference geodynamic solver. Use when the user wants to benchmark on Pyroclast Stokes & Advection Benchmarks, or asks about evaluating this task. Reports parallel speedup.
metadata:
  skill_kind: dataset_eval
  source_arxiv: 2603.14040
  bibtex_key: ferrari2026pyroclast
  confidence: high
---

# pyroclast-benchmark-eval

> Pyroclast: A Modular High-Performance Python Solver for Geodynamics — Ferrari (2026) (arXiv:2603.14040, 2026)

## What this evaluates

Evaluates the numerical correctness, stability, and parallel performance (CPU/GPU strong scaling and distributed weak scaling) of a matrix-free finite difference geodynamic solver.

## Datasets

- **Pyroclast Stokes & Advection Benchmarks** — total ?; splits: stokes-solver (-1), distributed-advection (-1); repo https://github.com/MarcelFerrari/Pyroclast

## Metrics

- `parallel speedup` **(primary)** — range: other
  - Defined as $S_p = T_1 / T_N$, where $T_1$ is the baseline runtime with a single thread/rank and $T_N$ is the runtime with $N$ threads/ranks. Measures scaling efficiency across hardware configurations.
- `relative energy residual` — range: [0, 1]
  - The solver's surrogate residual metric after 500 Uzawa iterations, indicating convergence to the physical solution. Target convergence is ~10^{-4}.
- `total runtime` — range: other
  - Summed execution time over multiple solver iterations after a warm-up phase, recorded in seconds. Used as the baseline for speedup calculations.

## Input / output format

**Input**: Structured grid configurations (2500x2500 to 15000x15000), material parameters, boundary conditions, and solver hyperparameters (6 grid levels, Jacobi smoother, 5 pre/post iterations, 500 Uzawa cycles). For distributed tests: MPI rank decomposition, local grid size per node (5000x5000), and marker advection settings.

**Output**: Numerical outputs including convergence trajectories (residual per iteration), total execution time in seconds, and computed parallel speedup values.

## Scoring recipe

```python
def compute_speedup(T1, TN):
    return T1 / TN

def compute_residual(final_res, initial_res):
    return final_res / initial_res

# Protocol:
# 1. Execute warm-up phase to stabilize memory/cache state.
# 2. Run solver for N iterations, record total time -> T1 (1 thread) or TN (N threads).
# 3. Speedup = T1 / TN.
# 4. Track relative energy residual over 500 Uzawa cycles.
# 5. Report mean speedup and final residual across repeated runs.
```

## Common pitfalls

- CPU scaling is memory-bound and stagnates or degrades beyond 32-64 threads due to NUMA/cache limitations, not algorithmic inefficiency.
- GPU speedup is highly problem-size dependent; grids smaller than 2000x2000 underutilize the GPU and yield low speedup (~2.4x).
- Results require explicit process pinning, NUMA binding, and a dedicated warm-up phase to avoid transient effects and ensure reproducibility.

## Evidence (verbatim from paper)

> Figure [20] reports the total runtime of 500 Uzawa iterations for each configuration, while Figure [20] shows the corresponding parallel speedup defined as $S_{p}\=T_{1}/T_{N}$, where $T_{1}$ is the runtime with 1 thread and $T_{N}$ with $N$ threads. Each benchmark included a dedicated warm-up phase preceding the timed execution. During the benchmark run, multiple iterations of the solver were performed, and the total runtime over all iterations was recorded.

## Citation

```bibtex
@misc{ferrari2026pyroclast,
  title={Pyroclast: A Modular High-Performance Python Solver for Geodynamics},
  author={Ferrari (2026)},
  year={2026},
  note={arXiv:2603.14040}
}
```

- arXiv: 2603.14040