Mixed-precision GPU-native reservoir simulation

01 · Overview

A fully implicit simulator that lives on the GPU

NVRS-CUDA is a three-phase black-oil simulator with dissolved gas and vaporised oil in which property evaluation, analytic Jacobian assembly, algebraic multigrid and the well equations all run on the device. Each Newton step is solved with flexible GMRES and a two-stage constrained-pressure-residual (CPR) preconditioner: stage 1 is one AmgX algebraic-multigrid V-cycle on the pressure block, stage 2 is 2 x 2 block Jacobi. The state, residuals and Krylov recurrences stay in double precision; the preconditioner runs in lower precision. BoomerAMG from HYPRE and Ginkgo can replace AmgX in stage 1.

GPU-resident Newton loop

The host reads the deck and writes results; inside a Newton iteration nothing returns to the host.

Mixed-precision CPR

Stage 1: one AmgX V-cycle on the pressure block in FP32. Stage 2: 2 x 2 block Jacobi. Outer FGMRES in FP64 keeps the answer exact.

Block-sparse operator, CUDA graphs

A block-sparse operator cuts the cost per Krylov iteration; graph replay is used only where it pays for itself.

Validated against OPM Flow

Norne, SPE10 model 2 and synthetic models of up to five million cells, well by well and field by field.

02 · Mathematics

Why reduced precision is safe here

In plain words

Newton handles the nonlinear conservation equations: each time step becomes a sequence of linear systems.
FGMRES solves each linear system, in FP64, and measures the true residual at every iteration.
CPR is the preconditioner that makes FGMRES fast: a multigrid V-cycle for the pressure part, then a cheap 2 x 2 block correction for everything else.
Reduced precision is used only inside the preconditioner, where it moves fewer bytes; FP64 FGMRES and the Newton tolerances decide when the solution is accurate enough.

Two-stage CPR: one AmgX V-cycle approximates the inverse of the pressure block A_p, then 2 x 2 block Jacobi (M_2) acts on the full system.

Flexible GMRES minimises the true residual over whatever directions the preconditioner produced, so rounding inside the preconditioner changes the iteration count, not the solution.

The condition that licenses an FP32 preconditioner: FP32 unit roundoff times the conditioning of the preconditioned operator.

03 · The simulator, live

Inside one Newton step: FGMRES with two-stage CPR

This is the simulator in miniature, running in your browser. A fully implicit two-phase model has pressure and water saturation in every cell, so the Jacobian has 2 x 2 blocks. Every time step is a Newton loop, and every Newton step solves J δu = −R with FGMRES in FP64, preconditioned by two-stage CPR. Stage 1 applies quasi-IMPES weights to form the pressure matrix and one aggregation-multigrid V-cycle (the role AmgX plays). Stage 2 is 2 x 2 block Jacobi. Choose the precision of the preconditioner and watch the iteration counts, the time per stage and the solution.

1Assemble R, Janalytic Jacobian, 2 x 2 blocksFP64

2FGMRESouter Krylov on the full systemFP64

3CPR stage 1quasi-IMPES pressure matrix, one multigrid V-cycleFP32

4CPR stage 22 x 2 block Jacobi on the full systemFP32

5Newton updatechop, converge, grow or cut the stepFP64

Preconditioner precision Wells Grid Heterogeneity (std of ln K) Oil/water viscosity ratio Geology seed FGMRES restart Time per movie frame End time Show

Interactive CPU implementation of the algorithm, in JavaScript in your browser. The GPU results on this page come from the separate research implementation.

Two-phase oil and water, Corey relative permeability, slightly compressible rock. In the five-spot, water is injected at a fixed rate in the centre and four producers at fixed bottom-hole pressure sit in the corners, with a Peaceman-type well index proportional to the well-cell permeability. Every run is repeated with an FP64 preconditioner as the reference (red); the curves of the chosen precision (blue) lie on it. The model has no gas phase, so there is no WGPR. Newton, the update limits and the time-step control use the same settings as the simulator in the paper: OPM Flow's convergence test (CNV below 1e-2 and material balance below 1e-6 in every phase), at most 20 Newton iterations, saturation changes limited to 0.2 per iteration, FGMRES(30) to a relative 1e-6, and dyadic chopping of each report step (up to six halvings, one released after two clean report steps). FP16 is emulated by rounding every operation of the preconditioner, applied to the row-scaled pressure matrix. This app runs in your browser, in JavaScript on the CPU, with no GPU and no server: a two-phase model far smaller than the three-phase GPU simulator, but with the same stages in its nonlinear and linear solvers. The run and its FP64 reference each run in a Web Worker, so the page stays responsive, and play as a movie: every frame advances the chosen time and redraws the field and the well curves. A fully implicit step costs far more per cell than the IMPES waterflood further down this page, so the 56 x 56 grid plays at about one frame per half second.

04 · Movies

Fields against OPM Flow

**Norne: NVRS-CUDA against OPM Flow**Pressure, water, oil and gas saturation over the field history: OPM Flow, NVRS-CUDA and their difference.

**SPE10 model 2: the flood**1.12 million cells over 1,000 days: pressure and saturations from OPM Flow and NVRS-CUDA, and their difference.

**Three-phase benchmark**A three-phase waterflood with free gas: OPM Flow, NVRS-CUDA and the difference at every report step.

05 · Results

Norne, SPE10 and the Million family

Wall clock in seconds on one GPU, from the paper. Flow is OPM Flow 2024.10 (its GPU linear solver on one rank for Norne and Bench3P, its CPU solver with two threads for SPE10 on the VR72 host). NVRS-CUDA runs every model on one GPU with BoomerAMG or AmgX in CPR stage 1.

Model	Cells	GB200			VR72
Model	Cells	OPM Flow	BoomerAMG	AmgX	OPM Flow	BoomerAMG	AmgX
Norne field	44,431	488.51	71.51	44.95	243.88	73.54	35.73
Bench3P	4,800	9.25	3.00	0.89	5.32	3.44	0.95
SPE10 model 2	1,122,000	n/r	99.22	53.06	766.16	91.74	41.21
Million 1M	1,000,000	n/r	26.23	21.17	n/r	19.94	14.87
Million 2M	2,000,000	n/r	49.94	39.38	n/r	38.23	27.17
Million 4M	4,000,000	n/r	81.87	77.67	n/r	60.97	52.33
Million 5M	5,000,000	n/r	102.32	95.61	n/r	76.65	63.77

n/r: not run.

How the headline speed-ups were measured

Claim	Model	NVRS-CUDA	Baseline	What is timed	Caveat
3.5x	Norne, 44,431 cells, full history	46.75 s, one GB200 GPU, AmgX in CPR stage 1	165.62 s, OPM Flow 2024.10, CPU CPR solver on 16 ranks, the fastest Flow configuration tested on the GB200 node	wall clock of the whole run, median of three	GPU against CPU ranks on the same node
18.6x	SPE10 model 2, 1,122,000 cells, 1,000 days	41.21 s, one VR72 GPU, AmgX in CPR stage 1	766.16 s, OPM Flow 2024.10, CPU solver on two threads on the VR72 host	wall clock of the whole run	a GPU against two CPU threads, not a hardware-equivalent comparison

Both runs reproduce OPM Flow's fields and well rates (see Accuracy below, with the per-well R2 on Norne).

The Million family on one, two and four GPUs

Strong scaling of synthetic three-phase models from one to five million cells. Switch machines, or show the speed-up of each back end against its own one-GPU time; hover a bar for its value.

AmgX, 1 / 2 / 4 GPUs (light to dark)BoomerAMG, 1 / 2 / 4 GPUs

Hover a bar for its value.

Cells	Stage 1	GB200				VR72
Cells	Stage 1	1 GPU	2 GPUs	4 GPUs	1 to 4	1 GPU	2 GPUs	4 GPUs	1 to 4
1 M	BoomerAMG	26.23	33.85	43.27	0.61x	19.94	28.88	49.61	0.40x
	AmgX	21.17	18.96	18.26	1.16x	14.87	14.02	16.17	0.92x
2 M	BoomerAMG	49.94	52.84	65.32	0.76x	38.23	43.83	72.11	0.53x
	AmgX	39.38	31.05	30.84	1.28x	27.17	21.24	23.62	1.16x
4 M	BoomerAMG	81.87	80.00	86.67	0.94x	60.97	60.96	85.21	0.72x
	AmgX	77.67	58.01	51.78	1.50x	52.33	36.58	41.14	1.27x
5 M	BoomerAMG	102.32	96.16	96.64	1.06x	76.65	77.49	93.30	0.81x
	AmgX	95.61	71.14	61.71	1.55x	63.77	43.91	42.26	1.51x

Reading the scaling

More GPUs do not always help. With BoomerAMG in stage 1 the 1 M and 2 M models get slower from one to four GB200 GPUs (0.61x and 0.76x) and only break even at 5 M (1.06x), while AmgX gains from 1.16x to 1.55x. Splitting a fixed problem across more GPUs shrinks the work per GPU, but the halo exchanges and the setup of the multigrid hierarchy do not shrink with it; small problems hit that strong-scaling limit first. It is a property of the problem size and the communication, not a defect of either library, and it is reported here as measured.

06 · Performance

Speed on GB200 and VR72

**Against OPM Flow**Wall clock on the three-phase benchmark, one GB200 each, with the time split into linear solve and the rest.

**Every model, one GPU**Wall clock per model and back end on GB200 and VR72, from 4,800 to 5 million cells.

**Scaling on GB200**One to four GPUs, AmgX and BoomerAMG, 4,800 to 5 million cells.

**Scaling on VR72**The same ladder on VR72.

**Reusing the multigrid hierarchy**Freezing the AMG hierarchy across Newton iterations cuts the setup time by 8.7x on Norne.

**Inside the linear solve**Where the time goes inside CPR on Norne: AMG setup and V-cycles, stage 2, the operator and orthogonalisation.

**Where the work goes**Wall clock by phase for OPM Flow and NVRS-CUDA on the benchmark.

07 · Precision

Mixed precision and emulated FP64

**Mixed-precision Newton-FGMRES-CPR**The pipeline of one Newton step and the precision of every stage.

**Where each precision is spent**Share of time in native FP64, emulated FP64 and FP32 on Norne.

**Roofline**Reservoir SpMV and AMG sit far into the memory-bound region: moving fewer bytes matters more than faster arithmetic.

**Measured convergence**FP64 against mixed precision and mixed precision with emulated FP64, on GB200.

**Convergence across grid sizes**FP64 (solid) against mixed precision (dashed): the curves coincide.

**Emulated FP64 accuracy**Double-word (df64) arithmetic from pairs of FP32 values carries about 48 significant bits against 53 for FP64 (unit roundoff about 2^-48), as bounded in the paper.

**Emulated FP64 in the solver**Emulated-FP64 mixed precision tracks true FP64; plain FP32 stalls.

**Path to further speed-up**Projected solve speed-up as the preconditioner's share of the per-iteration byte traffic grows.

08 · Accuracy

Norne and SPE10 against OPM Flow

**Norne permeability**The field model with its wells.

**SPE10 model 2 permeability**Six orders of magnitude of permeability on a log scale.

**Norne, final state**Pressure and saturations: OPM Flow, NVRS-CUDA and the difference.

**Norne over time**Fields at selected report steps.

**SPE10 over time**The flood through the fluvial channels.

**Oil rate of every producer**NVRS-CUDA against OPM Flow over the full history.

**Agreement per well**R2 of every producer's rates and pressures against OPM Flow.

09 · Precision lab

FGMRES with a reduced-precision multigrid preconditioner

The argument of the paper, run live. A heterogeneous pressure equation, the system that CPR's first stage solves, is handed to flexible GMRES, preconditioned by one aggregation-multigrid V-cycle, the role AmgX plays in stage 1 of CPR. With the outer FGMRES in FP64, running the V-cycle in FP32 or even FP16 changes how many iterations are needed, not the accuracy reached. When the outer recurrences drop to FP32 as well, the residual stalls near FP32 roundoff. Every curve is the true residual, evaluated in FP64.

FP64 FGMRES + FP64 V-cycle FP64 FGMRES + FP32 V-cycle (the paper's design) FP64 FGMRES + FP16 V-cycle FP32 FGMRES + FP32 V-cycle

Grid Permeability contrast FGMRES restart Maximum iterations