Research code for the MPhil thesis Optimization of Large Systems of Equations for GPUs. It investigates how high-resolution shock-capturing finite-volume solvers for the compressible Euler and ideal MHD equations map to multicore CPUs and a single NVIDIA GPU.
The codebase provides a C++/OpenMP reference implementation and a CUDA implementation, together with YAML configurations and Slurm scripts used for validation, benchmarking, Intel VTune profiling, and NVIDIA Nsight Compute profiling. The main research focus is the intrinsic performance of solver kernels on one GPU: data layout, coalesced access, shared-memory organization, register pressure, occupancy, and memory-latency hiding.
Research-code status — Solver and optimization choices are compile-time source settings, not YAML options. The checked-in CPU and GPU programs currently build different profiles, described below. Several Slurm scripts and output paths preserve the original CSD3 experiment environment and must be adapted before reuse.
- Current Checked-in Profiles
- Numerical Model
- CPU and GPU Implementations
- Repository Layout
- Requirements
- Build
- Quick Start
- Thesis Results
- Validation and Known Limitations
- HPC and Profiling Scripts
- Thesis Citation
- Third-party Software and Licensing
| Executable | Equations | Numerical Flux | Precision | Parallel Profile |
|---|---|---|---|---|
CPU |
9-variable ideal MHD with GLM divergence cleaning | HLL | double |
OpenMP + Vector Class Library SIMD |
GPU_FORCE_SharedFlux_Euler |
4-variable Euler | HLLC | double |
CUDA, 16x16 blocks, 2 resident blocks/SM |
The CUDA target name is historical. With the current constants it uses HLLC, not FORCE, and
USE_SHARED_FLUXisfalse. Always inspect the compile-time profile before interpreting an executable name or benchmark result.
Active settings are located in:
| File | Controls |
|---|---|
CPU.cpp |
CPU equation/state aliases and the active flux path |
GPU.cu |
GPU equation/state aliases |
includes/GPU/Constants.cuh |
Precision, Riemann solver, slope limiter, block geometry, launch bounds, fast division, profiling, and shared-flux switches |
Changing any of these settings requires rebuilding the corresponding target. The RiemannSolver and LimiterType keys at the top of the Shock Bubble YAML file are not consumed by the current executable.
Both implementations use a two-dimensional finite-volume method on a Cartesian grid. Each time step applies dimensional splitting in the x and y directions:
- Reconstruct piecewise-linear interface states with a slope limiter.
- Evolve the reconstructed states by a half time step.
- Compute interface fluxes with HLL, HLLC, or the SLIC/FORCE reference path.
- Apply a Godunov finite-volume update.
Available physics and numerical components
- Compressible Euler equations with an ideal-gas equation of state.
- Ideal MHD with a 9-component state:
[rho, rho*vx, rho*vy, rho*vz, Bx, By, Bz, E, Phi]. - Dedner-style GLM hyperbolic-parabolic cleaning for the MHD divergence constraint.
- HLL and HLLC approximate Riemann solvers.
- SLIC with the FORCE centered flux as a Riemann-solver-free reference.
- Minbee, Van Albada, Van Leer, and Superbee limiter implementations on the GPU.
HLLD-related source remains experimental and is not a supported solver in this snapshot.
| CPU | CUDA | |
|---|---|---|
| Data layout | Row-major C++ containers, two ghost-cell layers | Structure-of-arrays device layout for coalesced access |
| Parallelism | OpenMP over outer grid loops | One CUDA thread per cell within overlapping tiles |
| Vectorization | Vector Class Library SIMD in reconstruction and updates | — |
| Memory | Cached reconstructed states and fluxes | Pitched allocation (cudaMallocPitch), constant memory for simulation constants |
| Reduction | — | Thrust reduction for maximum wave speed |
| Kernel design | — | Fused x/y compute kernels, dynamic shared-memory tiles |
| Occupancy | — | __launch_bounds__ control of resident blocks per SM |
| Profiling | Intel VTune | NVIDIA Nsight Compute |
| Shared flux | — | Optional shared-flux kernels keeping the Godunov update on chip |
The CUDA implementation uses a stride-4 overlap between neighboring tiles so that each fused kernel can reconstruct interface states and update an interior region while retaining the stencil halo in shared memory.
.
├── CMakeLists.txt # C++/CUDA build and target definitions
├── CPU.cpp # OpenMP CPU driver
├── GPU.cu # CUDA driver
├── CONFIG_CPU/
│ ├── Euler/ # CPU Euler YAML and Slurm configurations
│ └── MHD/ # CPU MHD YAML and profiling scripts
├── CONFIG_GPU/
│ ├── Euler/ # GPU Euler YAML and Slurm configurations
│ └── MHD/ # GPU MHD YAML and Nsight scripts
└── includes/
├── CPU/ # CPU data structures and numerical methods
├── GPU/ # CUDA data structures, kernels, and methods
└── third_party/ # Vendored VCL, yaml-cpp, and inactive HighFive
- Linux or another environment supported by CUDA and OpenMP
- CMake ≥ 3.18 (the build file declares an older minimum, but it uses CUDA architecture handling introduced in 3.18)
- C++20 compiler with OpenMP support
- NVIDIA CUDA Toolkit with
nvccand Thrust - An NVIDIA GPU compatible with the selected CUDA architecture
The checked-in build fixes
CMAKE_CUDA_ARCHITECTURESand NVCC code generation to80/sm_80, matching the thesis-era NVIDIA A30 setup. Edit both settings inCMakeLists.txtfor a different GPU architecture.
Dependencies: Vendored yaml-cpp and the header-only Vector Class Library, plus OpenMP and CUDA/Thrust from the toolchain. HighFive and HDF5 are present as inactive research dependencies (CMake integration is commented out). Because the CMake project enables both CXX and CUDA, configuration requires a working CUDA compiler even when only the CPU target is requested.
From the repository root:
cmake -S . -B build
cmake --build build --parallelTo build a single target:
cmake --build build --parallel --target CPU
cmake --build build --parallel --target GPU_FORCE_SharedFlux_EulerThe following commands match the checked-in profiles and use the smallest provided configurations.
CPU — MHD / HLL:
OMP_NUM_THREADS=8 ./build/CPU CONFIG_CPU/MHD/OrszagTang.yamlGPU — Euler / HLLC:
./build/GPU_FORCE_SharedFlux_Euler CONFIG_GPU/Euler/CylindricalExplosion.yamlEach executable requires exactly one YAML path. The current drivers do not guard against a missing positional argument, so always invoke them with a configuration file.
YAML schema
| Key | Required by | Meaning |
|---|---|---|
Nx, Ny |
CPU and GPU | Interior grid resolution |
Xmin, Ymin |
Optional | Lower bounds; default to 0 |
Xmax, Ymax |
CPU and GPU | Upper bounds |
Tmax |
CPU and GPU | Final simulation time |
InitialCondition |
CPU and GPU | Named built-in test problem |
BoundaryCondition |
CPU and GPU | Transmissive, Periodic, or KelvinHelmholtz |
CFL |
CPU and GPU | Courant number |
Gamma |
CPU and GPU | Ratio of specific heats |
OutputTs |
CPU and GPU | Requested output times; currently inactive |
threshold |
CPU and GPU | Output-time tolerance; currently inactive |
OutputFolder |
CPU and GPU | Output path; currently inactive and often site-specific |
schedule.type |
CPU | Parsed OpenMP schedule name |
schedule.chunk_size |
CPU | OpenMP schedule chunk size |
The CPU parser accepts STATIC, DYNAMIC, GUIDED, and AUTO, but the current driver ultimately calls omp_set_schedule with guided. Treat the YAML schedule name as documentary until that hard-coded call is corrected.
Bundled test cases
| System | Cases |
|---|---|
| Euler | Cylindrical Explosion, Quadrant, Shock Bubble |
| Ideal MHD | Orszag-Tang vortex, MHD Rotor, Kelvin-Helmholtz |
| GPU-only MHD | Cloud Shock |
Profile YAML files under CONFIG_GPU/MHD/ use larger meshes intended for profiling runs, not quick functional checks.
Both drivers print progress and total wall-clock time to standard output. CSV snapshot blocks are commented out, so OutputTs, threshold, and OutputFolder are parsed but do not produce result files.
If CSV output is restored, update the thesis-era absolute OutputFolder paths in the YAML files first. The intended filenames are:
{OutputFolder}/{Nx}_{Ny}_{time}_CPU.csv
{OutputFolder}/{Nx}_{Ny}_{time}_GPU.csv
The thesis benchmarked multiple compile-time variants rather than only the defaults in this checkout. Headline speedups are peak values relative to the serial single-core CPU implementation, not relative to the OpenMP build.
| Optimized MHD Solver on One NVIDIA A100 | Peak Speedup vs. Single-Core CPU |
|---|---|
| HLLC | ~360x |
| HLL | ~433x |
| SLIC/FORCE | ~487x |
The shared-flux optimization improved the highest-resolution MHD runs by approximately:
| Solver | MHD Improvement | Euler Improvement |
|---|---|---|
| HLLC | ~30% | ~25–30% |
| HLL | ~33% | — |
| SLIC | ~35% | — |
The optimized path stores numerical fluxes in an additional shared memory array, reducing global-memory traffic during the Godunov update.
Additional findings
- A
16x16block was the best tested block geometry for these kernels. - Two resident blocks per SM gave the best MHD performance balance; requesting more caused costly register spilling.
- Pitched global-memory allocation added roughly 3–4% on large grids.
- Fast reciprocal-based division added roughly 5–10% in some runs but caused a high-resolution HLLC case to lose convergence, so it was excluded from the final optimized profile.
- Observed branch efficiency remained above 85%; within that measured range, branch divergence did not correlate with a material runtime change. Memory latency, register pressure, and shared/global-memory traffic dominated.
- The 32-thread OpenMP MHD implementation achieved less than 15x maximum speedup and became DRAM-bandwidth limited in the y-direction on large grids.
Benchmark context
| Hardware / Toolchain | |
|---|---|
| Euler CPU baseline | Two Intel Xeon Silver 4314, 32 physical cores |
| Euler GPU | One NVIDIA A30 (two-GPU node) |
| MHD GPU profiling & final speedups | NVIDIA A100 |
| Euler toolchain | GCC 13.3.0, NVCC 12.9.41 |
| Reporting | Averages of five runs (Euler) |
Euler and MHD results used different GPUs, so their relative speedups should not be attributed to equation complexity alone. Speedup against a serial baseline is also sensitive to the quality of that baseline; the thesis recommends adding hardware-independent metrics such as FLOP/s in future work.
The thesis validated the CUDA solver against the CPU implementation and published test-case behavior:
- Euler: CPU/GPU normalized L2 differences remained below
1e-10in the reported tests. - Ideal MHD: Differences were larger (~
1e-4to1e-2), and some high-resolution Orszag-Tang and Kelvin-Helmholtz variants diverged or produced NaNs. The thesis attributes this partly to floating-point evaluation order, fused operations, and MHD sensitivity, but does not claim the issue is fully resolved.
Known limitations of this research snapshot
- No project-level automated test suite or CI workflow.
- The default CUDA profile does not enable the final shared-flux optimization.
- CPU Euler switching is experimental: part of the generic reconstruction path still contains MHD-specific access to the
Phicomponent. - The CPU/GPU time step is computed as
CFL * Dx / amax; it does not usemin(Dx, Dy)or clamp the last step toTmax/an output time. - HLLD was not completed.
- The two-dimensional directional split is first order in its time composition.
- Fast division remains disabled because of the observed stability trade-off.
- The current framework is single-material and has no adaptive mesh refinement.
The .sbatch files under CONFIG_CPU/ and CONFIG_GPU/ document the original CSD3 experiment workflow. They are archival reproduction aids rather than portable launch scripts.
Before submitting, review:
- Slurm account, partition, module, time, CPU, and GPU directives.
- Hard-coded
/home/kw623/.../GPU_Optimizationpaths. - Historical executable names such as
CPU_HLLC,GPU, orGPU_FORCEthat are not produced by the current CMake file. - Output and temporary-file locations.
Intel VTune was used for CPU profiling and NVIDIA Nsight Compute for CUDA kernel profiling. Generated CSV, log, image, and .ncu-rep artifacts are ignored by the top-level Git rules.
If this repository supports academic work, please cite:
Keyuan Wang, Optimization of Large Systems of Equations for GPUs, MPhil dissertation, Cavendish Laboratory, University of Cambridge, February 2026.
Vendored dependencies retain their own license files under includes/third_party/. This repository currently has no top-level project license; do not assume permission beyond the licenses of those dependencies.