# Porting of ISSM to GPUs

**URL:** <https://forum.access-hive.org.au/t/porting-of-issm-to-gpus/6656>\
**Category:** ACCESS-ISSM\
**Created:** [17 July 2026 03:51 UTC](https://forum.access-hive.org.au/t/porting-of-issm-to-gpus/6656 "2026-07-17T03:51:40Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![justinh2002](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/justinh2002/32/2853_2.png) [@justinh2002](https://forum.access-hive.org.au/u/justinh2002)\
**Post date:** [17 July 2026 03:51 UTC](https://forum.access-hive.org.au/t/porting-of-issm-to-gpus/6656/1 "2026-07-17T03:51:41Z")

</div>

# ISSM GPU Solver — Change Notes

Work to make ISSM run its linear solves on NVIDIA V100 GPUs (NCI Gadi  
`gpuvolta` queue) via a CUDA-enabled PETSc. Chronological summary of what was  
changed and why.

## Files changed / added

| File | Change |
| --- | --- |
| `externalpackages/petsc/install-3.22-gadi-gpu.sh` | Build PETSc against **system OpenMPI** instead of downloading MPICH |
| `configure.sh` (gpu step) | Added `--with-python=...` so Python wrappers build |
| `src/m/classes/gpuoptions.py` | New solver toolkit: CUDA vec/mat + GMRES/GAMG, with block-size and host-PtAP tuning |
| `src/c/Makefile.am` | Added missing `AndersonAccelerator.cpp` to the build |
| `test/gpu_kspsolve_test.py` | Resolution-sweep CPU-vs-GPU benchmark |
| `test/run_gpu_kspsolve_test.pbs` | PBS job script for the gpuvolta queue |
| `test/gpu_solver_report.tex` | LaTeX write-up of the work |

## Problems hit and how they were fixed (in order)

1. **`ModuleNotFoundError: IssmConfig_python`** — GPU `configure` was missing  
`--with-python`, so the compiled Python wrappers (`*.so`) were never built.  
Added `--with-python=/apps/intel-python3/.../python3`. Wrappers install to  
`$ISSM_DIR/lib`, which must be on `PYTHONPATH`.

2. **`libk5crypto.so.3: undefined symbol EVP_KDF_ctrl`** — Intel Python ships an  
old `libcrypto` lacking `EVP_KDF_ctrl`. Fixed at _runtime_ in the PBS script  
by `LD_PRELOAD`-ing system OpenSSL 1.1.1k (`/lib64/libcrypto.so.1.1`,  
`libssl.so.1.1`). Not a build change.

3. **Tiny problem, GPU slower than CPU** — expected. At 50 km (680 DOFs) MUMPS  
is sub-second; GPU overhead (CUDA init, CUSPARSE assembly, H2D transfer)  
dominates. Built the resolution-sweep benchmark to find the crossover.

4. **GMRES + bjacobi/ILU stalls at 5 km+** — ILU preconditioning quality decays  
with mesh refinement. Switched default `pc_type` from `bjacobi` to **`gamg`**  
(algebraic multigrid).

5. **`np=4` multi-rank → 4 independent rank-0 jobs** — root cause: `issm.exe`  
was linked against **two MPI runtimes at once** (`libmpi.so.12` MPICH from  
GPU PETSc + `libmpi.so.40` OpenMPI from `mpicxx`). Harmless at `np=1`, breaks  
any `np>1`. **Fixed by rebuilding GPU PETSc against system OpenMPI** : removed  
`--download-mpich=1`, added `--with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90` + an MPI sanity check in the install script. Requires  
`rm -rf install-gpu src-gpu` before re-running `source configure.sh gpu`.

6. **Link error `undefined reference to AndersonAccelerator::...`** — after the  
clean rebuild, `AndersonAccelerator.cpp` was missing from `src/c/Makefile.am`  
(the file existed and was referenced by `solutionsequence_nonlinear.cpp` but  
never compiled). Added it to the source list.

7. **GAMG won’t converge at 2 km / 1 km** — GAMG treated the SSA system as a  
scalar matrix. SSA has **2 coupled DOFs/node (vx, vy)**. Added  
`mat_block_size = 2` so GAMG aggregates by node. Fixed convergence for all  
resolutions that fit in memory.

8. **`np=4` on a single GPU is ~4× slower** — 4 ranks contend for one V100 and  
add MPI overhead inside GAMG with no parallelism gain. Reverted benchmark to  
`NP = 1`. (Multi-rank only helps with one GPU per rank; the MPI fix is still  
needed for that future case.)

9. **2 km / 1 km → `cudaErrorMemoryAllocation` in `MatProductSymbolic`** —  
GAMG’s coarse-grid Galerkin product (`MatPtAP`) uses CUSPARSE SpGEMM, whose  
scratch buffers overflow 32 GB at 397k / 1.59M DOFs. Fix is to run the matrix  
products on the host (96 GB) while the Krylov iteration stays on the GPU —  
but **the option name matters** : GAMG calls `MatPtAP()` directly  
(`api_user=true`), so PETSc reads `-matptap_backend_cpu`, _not_ the generic  
`-mat_product_algorithm_backend_cpu` (which it silently ignored —  
`Option left: ... value: true`). Set `matptap_backend_cpu = true` plus  
`matmatmult_backend_cpu = true` (for the inner A·B products of the decomposed  
PtAP). This cleared the OOM: **2 km then passed at 2.2×**.

10. **1 km → true residual stuck at 1.51e-6 (\> 1e-6 gate)** — _not_ OOM and _not_  
iteration-starved: the residual was **bit-identical** across `ksp_max_it`  
2000→5000 and `ksp_gmres_restart` 30→100, proving GMRES was converging to its  
relative tolerance, not exhausting iterations. `ksp_rtol = 1e-10` (relative  
reduction of the _preconditioned_ residual) left the _true_ residual  
`||KU-F||/||F||` at 1.5e-6. Tightened `ksp_rtol = 1e-12` → true residual  
~1.5e-8. **1 km then passed at 3.2×.** (Longer `ksp_gmres_restart = 100` kept  
so the extra iterations to reach the tighter tol don’t stagnate.)

## Multi-GPU (one rank per GPU) — the parallel saga

With single-GPU complete, the next goal was multi-GPU scaling (1 MPI rank per  
V100 on a 4-GPU gpuvolta node). PETSc maps rank `r` → GPU `r % ndev`  
automatically (`PETSC_DECIDE`, verified in `cupmdevice.cxx:259`), so no binding  
wrapper is needed. Getting it to actually converge took a four-step hunt:

1. **CPU MUMPS scaled, every GPU solve failed.** First confirmation the OpenMPI  
rebuild works at np\>1: 4-rank MUMPS sped up (1km: 865→395 s). GPU solves  
crashed in `PCGAMGOptimizeProlongator_AGG`: _“Computed maximum singular value  
as zero.”_ (smoothed-aggregation eigenvalue estimate).

2. **`pc_gamg_agg_nsmooths` silently ignored.** Setting it to integer `0` hit the  
`toolkits.py` marshaller bug (`if not optionvalue:` writes a _valueless flag_  
for falsy values), so PETSc kept its default. **Pass option values as  
STRINGS** (`'0'`) to survive marshalling. (Latent foot-gun for any option set  
to integer 0.)

3. **Unsmoothed aggregation then _diverged_** (residual 3–7, worse than the RHS).  
A `-use_gpu_aware_mpi 0` guess changed nothing (bit-identical residuals), so  
that was a red herring — ruling it out mattered.

4. **Isolation diagnostic pinned the real cause** (`gpu_mgpu_diagnose.py`: one  
small mesh, np=4, several PCs back-to-back in one job). Results:  
`bjacobi+ilu` CONVERGED (so the **parallel GPU matvec is correct** ); GAMG with  
**host-PtAP FAILED** but the _identical_ GAMG with **GPU-PtAP CONVERGED  
(1e-11) **. So** `matptap_backend_cpu=true` builds a corrupt coarse operator in  
parallel** — the very flag that fixed the single-GPU OOM. The rule:  
**host-PtAP for single-GPU (memory); GPU-PtAP for multi-GPU (correctness, and  
it fits since each rank holds ~1/N of the matrix).** Encoded as an NP-based  
switch in `gpu_kspsolve_test.py`.

### Multi-GPU benchmark (4× V100, unsmoothed aggregation + GPU-PtAP) — all PASS

| Resolution | DOFs | 4-GPU (s) | 1-GPU (s) | 1→4 GPU scaling |
| --- | --- | --- | --- | --- |
| 5 km | 64,362 | 9.58 | 10.42 | 1.09× |
| 2 km | 396,752 | 40.61 | 63.30 | 1.56× |
| 1 km | 1,585,078 | 149.85 | 289.49 | **1.93×** |

Scales sublinearly (~48% efficiency at 1.6M DOFs, growing with size — small  
problems are communication-bound). NB the 1→4 comparison mixes a preconditioner  
change (multi-GPU used unsmoothed, single-GPU smoothed); a clean strong-scaling  
study needs the same config at np=1,2,4.

## A100 benchmark and extended investigation

### 15. Three-column benchmark on A100 (July 2026)

Submitted the same three-solver sweep to the `dgxa100` queue (single A100 80GB,  
`ISSM_GPU_PTAP=1`). With 80GB HBM2e, the CUSPARSE SpGEMM scratch for  
`MatPtAP` fits without OOM (unlike 32GB V100), so GPU-PtAP can run on the device.

Result: GPU-PtAP converges on all resolutions up to 1km — **but only with  
 `pc_gamg_agg_nsmooths=0`** (unsmoothed aggregation). Smoothed aggregation  
(`nsmooths=1`) diverges: `solver residue too high: 5.38 > 1e-6`. Root cause  
identified later (see #20).

| Resolution | DOFs | MUMPS (s) | CPU-GMRES (s) | GPU-GMRES (s) | vs MUMPS | vs CPU-GMRES |
| --- | --- | --- | --- | --- | --- | --- |
| 50 km | 680 | 2.638 | 2.171 | 2.697 | 0.98× | 0.80× |
| 20 km | 4,172 | 2.525 | 2.296 | 3.350 | 0.75× | 0.69× |
| 10 km | 16,304 | 4.741 | 4.494 | 4.883 | 0.97× | 0.92× |
| 5 km | 64,362 | 16.055 | 12.815 | 12.397 | 1.30× | 1.03× |
| 2 km | 396,752 | 135.062 | 81.540 | 67.680 | 2.00× | 1.20× |
| 1 km | 1,585,078 | 821.372 | 374.803 | 267.222 | **3.07×** | **1.40×** |

Hardware speedup improved to **1.40×** (vs 1.35× on V100) because GPU-PtAP  
eliminates the H2D transfer tax that dominated on V100 (confirmed by nsys: only  
~2% of V100 walltime was on GPU with host-PtAP). PBS reported 16% GPU  
utilisation with GPU-PtAP vs ~2% with host-PtAP.

### 16. Solver variants sweep (Chebyshev, BiCGStab, pipelining)

Tested five options combinations on A100 NP=1 at 2km and 1km  
(`gpu_speedup_variants_test.py`):

| Variant | 1km (s) | vs MUMPS | vs CPU-GMRES |
| --- | --- | --- | --- |
| baseline (GMRES, SOR) | 300.4 | 2.88× | 1.37× |
| Chebyshev+Jacobi smoother | 303.6 | 2.84× | 1.35× |
| BiCGStab | 296.5 | 2.91× | 1.38× |
| BiCGStab + Chebyshev | **294.2** | **2.94×** | **1.40×** |
| nsmooths=1 + CG eigenestimator | FAIL | — | — |

All within ~6% of each other. No variant breaks through the ~1.40× hardware  
ceiling. The CG-eigenestimator fix did not rescue nsmooths=1 (root cause is  
different — see #20). **cheb+bcgs is marginally best at 1km** but the gain is  
too small to justify changing the default.

### 17. Multi-GPU A100 (NP=4)

Ran the benchmark with NP=4 (one rank per A100) on the dgxa100 node.  
`NP > 1` triggers GPU-PtAP + nsmooths=0 automatically in `gpu_kspsolve_test.py`.

| Config | 1km GPU-GMRES | 1km CPU-GMRES | vs CPU-GMRES | vs MUMPS |
| --- | --- | --- | --- | --- |
| A100 NP=1 | 267s | 375s | **1.40×** | 3.07× |
| A100 NP=4 | 182s | 177s | **0.97×** | 2.14× |

Going from NP=1 to NP=4: CPU-GMRES scales 2.1× (375→177s) but GPU-GMRES only  
scales 1.5× (267→182s). At NP=4 the GPU barely matches CPU. Reason: with 4  
ranks each holding ~400k DOFs, the GMRES global AllReduce (one per iteration)  
and halo exchanges dominate. GPU is ~64% utilised but spending most of its time  
waiting on MPI rather than computing. PBS: GPU utilisation 64%, GPU memory only  
6.35 GB (1/4 of the single-GPU problem).

**Single-GPU is the sweet spot** for this problem size. Multi-GPU is  
communication-limited, not compute-limited.

### 18. GPU-aware MPI (NP=4)

OpenMPI 4.1.3 on Gadi is compiled with `--with-cuda` (`smcuda` BTL,  
`MCA coll: cuda`), so GPU-aware MPI can be enabled by removing  
`-use_gpu_aware_mpi 0`. Tried this on NP=4 to eliminate the  
GPU→CPU→MPI→CPU→GPU round-trip per iteration.

Result: **GPU GMRES fails at 2km and 1km** (residual 42026 at 2km, 1.43 at 1km).  
Small resolutions pass (GPU-PtAP not triggered). The GPU-aware MPI transfers  
corrupt the parallel coarse operator assembly during GPU-PtAP. UCX’s CUDA IPC  
path cannot safely share device pointers across CUDA contexts on different GPUs.

**GPU-aware MPI is incompatible with GPU-PtAP at NP\>1.** Keep  
`-use_gpu_aware_mpi 0` for all multi-GPU runs.

### 19. Pipelined Krylov (NP=4, GPU-aware MPI off)

Pipelined variants (`pipefgmres`, `pipebcgs`) restructure the Krylov recurrence  
to overlap AllReduce(k) with SpMV(k+1), hiding the MPI latency.

| Variant | 2km (s) | 1km (s) | vs CPU-GMRES |
| --- | --- | --- | --- |
| gmres (baseline) | 52.2 | 183.8 | 2.29× |
| bcgs | 52.9 | 184.0 | 2.29× |
| pipefgmres | 56.6 | 193.0 | 2.18× |
| pipebcgs | 54.6 | 191.6 | 2.20× |

Pipelining is **slightly slower** than standard GMRES at NP=4. On a single DGX  
node the 4 A100s are connected by NVLink (~600 GB/s), so the AllReduce latency  
is already negligible. The pipeline overhead (extra synchronisation checkpoints  
in the recurrence) outweighs any latency hiding. Pipelining would only help for  
inter-node communication (high-latency Infiniband), not intra-node NVLink.

Note: the CPU-GMRES baseline in this test was inadvertently set to nsmooths=0  
(same as GPU), making it appear slower (421s at 1km). The apparent 2.29×  
hardware speedup is misleading — it measures GPU vs CPU at the same (weaker)  
preconditioner quality, not best-vs-best. The true hardware speedup at NP=4  
remains ~0.97× (GPU-nsmooths=0 vs CPU-nsmooths=1).

### 20. nsmooths=1 GPU-PtAP root cause (CUSPARSE\_STATUS\_INSUFFICIENT\_RESOURCES)

Definitive diagnosis via `gpu_nsmooths_diag.py` — five variants at A100 NP=1,  
sweeping GPU/CPU assignment for MatMatMult (prolongator smoothing: A×P₀) and  
MatPtAP (coarse operator: P^T×A×P) independently.

The PETSc error log:

```auto
CUSPARSE_STATUS_INSUFFICIENT_RESOURCES (error 11)
MatProductSymbolic_SeqAIJCUSPARSE_SeqAIJCUSPARSE() ← GPU SpGEMM
MatProductSymbolic_ABC_Basic()
MatPtAP()
PCGAMGCreateLevel_GAMG()

```

**Root cause:** with nsmooths=1, the smoothed prolongator P has more non-zeros  
than the tentative P₀ (nsmooths=0). The CUSPARSE SpGEMM symbolic phase for  
P^T×A×P then exceeds CUSPARSE’s internal hardware resource limit (not device  
memory — a workspace/occupancy limit in the SpGEMM kernel). This is a hard  
CUSPARSE constraint that cannot be worked around via command-line options.

Isolation results:

- `n1-all-gpu` (GPU-MatMul + GPU-PtAP, n=1): FAIL (CUSPARSE crash → residual 5.38)
- `n1-mul-cpu` (CPU-MatMul + GPU-PtAP, n=1): FAIL (same CUSPARSE crash in MatPtAP)
- `n1-ptap-cpu` (GPU-MatMul + CPU-PtAP, n=1): ~926s at 2km (extremely slow but passes)
- `n1-all-cpu` (CPU-MatMul + CPU-PtAP, n=1): **49s at 2km (PASS, only 7% slower than n0-gpu)**

The fix for n1 is to put **both** MatMatMul and MatPtAP on CPU. At 2km this adds  
only 7% overhead vs n0-gpu because the better preconditioner (nsmooths=1) needs  
fewer GMRES iterations. This is essentially the same trade-off as the original  
V100 host-PtAP approach: CPU setup overhead offset by fewer Krylov iterations.

**The GPU-PtAP improvement and the nsmooths degradation cancel each other out:**  
both V100-host-PtAP+nsmooths=1 and A100-GPU-PtAP+nsmooths=0 land near 3.07–3.13×  
vs MUMPS at 1km. This is a fundamental result: the CUSPARSE SpGEMM workspace  
limit prevents smoothed aggregation from running on GPU with large matrices,  
and the preconditioner quality loss from nsmooths=0 offsets the GPU-PtAP gain.

## Final `gpuoptions` defaults

```auto
toolkit = petsc
vec_type = cuda
mat_type = aijcusparse
ksp_type = gmres
pc_type = gamg
mat_block_size = 2 # SSA: 2 DOFs/node (GAMG aggregation)
matptap_backend_cpu = true # GAMG PtAP on host (avoid GPU OOM); api_user name
matmatmult_backend_cpu = true # inner A*B of decomposed PtAP, also on host
mat_product_algorithm_backend_cpu = true # belt-and-braces (MatProduct path)
pc_gamg_coarse_eq_limit = 2000
pc_gamg_threshold = 0.08
ksp_rtol = 1e-12 # tight enough that TRUE residual clears 1e-6 gate
ksp_max_it = 5000
ksp_gmres_restart = 100

```

These defaults target **single-GPU**. For **multi-GPU** the driver overrides PtAP  
to run on the GPU (host-PtAP is broken in parallel):  
`matptap_backend_cpu=false`, `matmatmult_backend_cpu=false`,  
`mat_product_algorithm_backend_cpu=false`. Pass all option values as **strings**  
(toolkits.py marshalling drops falsy/zero values written as bare flags).

## Benchmark status (single V100, SquareShelf SSA, np=1) — COMPLETE, all PASS

Three-solver sweep isolates the hardware and algorithm contributions separately:

| Resolution | DOFs | CPU MUMPS | CPU GMRES+GAMG | GPU GMRES+GAMG | vs MUMPS | vs CPU-GMRES |
| --- | --- | --- | --- | --- | --- | --- |
| 50 km | 680 | 1.44 s | 0.84 s | 2.25 s | 0.64× | 0.37× |
| 20 km | 4,172 | 1.14 s | 1.21 s | 1.85 s | 0.62× | 0.65× |
| 10 km | 16,304 | 3.24 s | 3.04 s | 3.16 s | 1.03× | 0.96× |
| 5 km | 64,362 | 14.14 s | 12.54 s | 9.99 s | 1.42× | 1.25× |
| 2 km | 396,752 | 132.36 s | 78.16 s | 59.87 s | 2.21× | 1.31× |
| 1 km | 1,585,078 | 854.94 s | 368.82 s | 272.73 s | **3.13×** | **1.35×** |

The CPU GMRES+GAMG column uses identical GAMG tuning (`mat_block_size=2`,  
`ksp_rtol=1e-12`, etc.) with only `vec_type=standard` and `mat_type=aij` (no  
CUDA) — so the two rightmost columns change exactly one variable: the device.

**Reading the numbers:**

- **Algorithm effect (MUMPS → GMRES+GAMG on CPU):** ~2.3× at 1.6M DOFs. The  
iterative O(N) solver beats direct O(N^{3/2}) on the CPU alone. This is the  
dominant effect and grows with N.
- **Hardware effect (CPU → GPU, same GMRES+GAMG):** ~1.35× at 1.6M DOFs. The  
true GPU contribution — and it will grow further with the next resolution level  
as the bandwidth advantage of HBM2 over DRAM compounds.
- **Combined (MUMPS → GPU-GMRES, what a user experiences):** 3.13×.

Correctness: GPU vs CPU-MUMPS `Vx` relative diff ≤ ~4e-11 throughout.

## Rebuild / run cheatsheet

```bash
# interactive GPU node
qsub -I -q gpuvolta -l ncpus=12,ngpus=1,mem=96GB,walltime=2:00:00,storage=gdata/au88+scratch/ek9
module load cuda/12.3.2
export ISSM_DIR=/g/data1b/au88/jh7060/ISSM
cd $ISSM_DIR

# only needed if PETSc itself must be rebuilt (e.g. MPI/CUDA change):
rm -rf externalpackages/petsc/install-gpu externalpackages/petsc/src-gpu
source configure.sh gpu

# sanity check: must show ONLY libmpi.so.40 (OpenMPI), no libmpi.so.12
ldd $ISSM_DIR/bin/issm.exe | grep -i mpi

# run the benchmark (from a login node)
qsub $ISSM_DIR/test/run_gpu_kspsolve_test.pbs

```

Note: edits to `gpuoptions.py` or the test script are **Python-only** — no  
rebuild needed, just resubmit the PBS job. Only C/PETSc/Makefile.am changes  
require `source configure.sh gpu`.