Add a CPRW pressure stage to the system solver by hnil · Pull Request #7278 · OPM/opm-simulators

hnil · 2026-08-04T12:26:45Z

Supersedes #7209 (agreed with @ElyesAhmed), which implements the same feature by
bridging the system blocks into the classic PressureBhpTransferPolicy. This takes
the other route: assemble the coarse system directly from the B/C/D blocks
addBCDMatrix already extracts, so the pressure stage never reaches back into
well-model machinery.

The pressure stage of --linear-solver=system_cpr was reservoir-only, so the well
unknowns never reached the coarse system. This adds a CPRW stage for it, selected
with --linear-solver=system_cprw. The well DOF layout and the well weights are
produced in ISTLSolverSystem; below that point the stage sees only sparse
matrices, weights and integers. It defaults to the standard solver's weighting —
trueimpes on the reservoir equations, perforated-cell average on the well
equations.

Full NORNE_ATW2013, linear iterations and wall time:

np	`cprw`	`system_cpr`	`system_cprw`
1	2262 / 138s	4303 / 166s	2248 / 137s

With --convert-to-multisegment-well=per-connection:

	`cprw`	`system_cpr`	`system_cprw`
linear iterations	2716	4839	2558

Results agree with cprw to ~1e-5 on small cases and ~0.1% on Norne (different
Newton paths), serial vs parallel to ~5e-5. Norne also exercises well
opening/closing, i.e. the coarse-system rebuild.

preconditioner.well_transfer selects how the well unknowns take part in the
transfer: classic (default) is the classic cprw formulation, no_prolongation
restricts the well residual and leaves the segment correction to the trailing well
solve, full also prolongs the coarse bhp correction. A settings sweep on full
Norne picked classic, but by 2558 against 2569 -- the other two are there for
the inexact-well-solve case, where restricting the well residual should matter more
than it does with the exact solve used here.

Parity with #7209. Same: system_cprw option, extended (nCells + nWells)
pressure system, quasi-IMPES and cell-average well weights, no-perforations-on-rank
regularisation, zero-diagonal fallback, sequential and parallel, rebuild on
well-structure change, JSON validation and fixtures. Deliberately different: one
coarse unknown per well rather than per D block row, and the well weights are
handed down from the outer layer instead of being derived inside the
preconditioner. Config convention also differs — preconditioner.type stays
system_cpr with reservoir_solver.preconditioner.add_wells = true, mirroring how
cpr/cprw relate.

Still draft: the general system preconditioner (a TwoLevelMethodCpr and transfer
policy on the system operator, and a wrapper turning reservoir preconditioner + well
solver into a system preconditioner) is not in this PR.

🤖 Generated with Claude Code

ElyesAhmed

The SystemCprwPressureStage cleanly separates the well model from the linear‑solver core — WellDofLayout is simply a prefix‑sum data structure populated by ISTLSolverSystem, and the stage reuses the classic CPRW coarse‑operator types together with extendCommunicatorWithWells. As a result, the coarse solver and the parallel DOF convention remain consistent with the existing path.

The formulation is Ok and well documented, and the test coverage is fine too, particularly the independent dense R·S·P reference check. I also verified the singular/degenerate guards and the setupPropertyTree validation logic (requiring coarsesolver when add_wells is active, and mandating a flexible outer solver when a Krylov well solver is used), both of which are a fine.

I left a couple of minor inline notes — the well_transfer fallback being "full" vs. the shipped "classic", and the coarse‑matrix dump keying off preconditioner.verbosity which the generated config never sets. These are papercuts, not blockers.
Good to merge, but before #7289.

hnil · 2026-09-17T14:15:07Z

jenkins build this please

hnil · 2026-09-18T07:25:53Z

jenkins build this please

ElyesAhmed · 2026-09-18T11:32:43Z

jenkins build this please

ElyesAhmed

After this round, I will approve it!

The system solver (--linear-solver=system_cpr) solves the coupled (reservoir, well) system [A C; B D] explicitly, but its pressure stage was still reservoir-only: stage 1 of SystemPreconditioner ran a plain CPR on A with add_wells = false, so the well unknowns never entered the coarse system. That is exactly the coupling CPRW exists to supply. Classic CPRW cannot be reused here. PressureBhpTransferPolicy calls back into the live well model from inside the preconditioner (addWellPressureEquations -> per-well extractCPRPressureMatrix, which needs WellState just to ask isPressureControlled). That is what ties WellModelAsLinearOperator to TypeTag and what keeps CPRW out of NLDD and the GPU path. Instead the coarse system is assembled from the B/C/D blocks the outer layer already extracts, once, via addBCDMatrix. Two pieces of plain data are added to carry what is still missing: - WellDofLayout: a prefix sum over the per-well D dimensions, so the merged well block rows can be mapped back to wells (one row per standard well, one per segment for a multisegment well), plus the index of the pressure-like unknown inside a well block. - the well half of the weights. The weights calculator was already std::function<SystemVector()>; it now fills w[_1] as well as w[_0]. ISTLSolverSystem computes those weights (quasi-IMPES from the diagonal D block by default, 'unit' for debugging) and can later obtain them from the well model without the core changing. Both are produced in ISTLSolverSystem. Below that point the pressure stage reads nothing but sparse matrices, weights and integers -- no part of the well model is visible. The preconditioner factory previously discarded w[_1]; it now passes the whole SystemVector through and SystemPreconditioner derives the reservoir-only calculator for its sub-solvers. SystemCprwPressureStage builds the (Nres + nWells) scalar system as R*S*P with R = blockdiag(w0^T ; sum over the well's block rows of w1^T) and P = blockdiag(e_p ; e_q on the well's top block row), reusing Details::CoarseOperatorType and extendCommunicatorWithWells from PressureBhpTransferPolicy. One coarse unknown per well, as in classic CPRW; one per segment is left as a separate question. Two details differ from the classic policy: the well residual is really restricted, and the coarse well correction is prolonged back rather than discarded. Both are selectable through preconditioner.well_transfer, whose 'classic' value reproduces the classic formulation so that the two differ only in numerics. On SPE1 that lands on the classic cprw iteration count exactly. Selected with --linear-solver=system_cprw, or by setting preconditioner.reservoir_solver.preconditioner.add_wells in a JSON configuration. Linear iterations, serial: deck cprw system_cpr system_cprw SPE1 443 571 429 SPE9_CP 439 682 470 BASE2_MSW_HFA 30 82 52 and on SPE9_CP_SHORT the count stays flat under decomposition (110/107/105/106 for 1/2/4/8 ranks) where system_cpr drifts up (177/178/179/180). Also reject an approximate (Krylov) well solver combined with a non-flexible outer solver: such a well solve stops on a tolerance, so the preconditioner varies between applications and bicgstab or plain gmres are no longer valid. flexgmres is required. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit 6dd1454)

Two pieces of OPM#7209 that were missing here. The classic use_well_weights = false weighting is now available as well_weight_type = cellavg: average the reservoir weights over the cells a well block perforates and apply them to the conservation equations only, with weight zero on the control equation, falling back to unit weights when a well has no perforations on this rank. It is not a strictly worse or better choice than the quasi-IMPES default, so both are kept -- on SPE9_CP it needs 428 linear iterations against 470 for quasi-IMPES (and 439 for the classic cprw), while on SPE1 quasi-IMPES wins with 429 against 455. A well whose contraction cancels exactly leaves a zero on the coarse diagonal and makes the pressure system singular. Regularise that row to a unit diagonal instead of handing AMG a singular system. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit bf21d5d)

The pressure stage now defaults to well_weight_type = cellavg, so together with the trueimpes reservoir weights it already used, its weighting matches what cprw does by default (use_well_weights = false). Quasi-IMPES well weights remain available and are better on some cases, but the sensible default is the one the standard solver uses. A test pins both halves against setupCPRW so they cannot drift apart. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit dce7001)

The restriction and prolongation grow caller-supplied-output forms, so that a caller which already owns the target vectors does not need a copy, and the two existing forms delegate to them instead of repeating the loops. The coarse matrix, communication and rhs also get handles, so the coarse level can be driven from outside the stage. No functional change: the transfers compute what they computed before, and the stage is still the only thing that assembles the coarse system. (cherry picked from commit 66d642a, with the composable preconditioner parts and the Dune transfer policy left out -- they belong with the general system preconditioner, which is not part of this PR.) Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

For a multisegment well the coarse column was taken from the well's top block alone, while MultisegmentWellEquations::extractCPRPressureMatrix accumulates over every segment row. With one segment per connection that discarded most of the well: on Norne with --convert-to-multisegment-well=per-connection the coupled pressure stage needed 4939 linear iterations against 2716 for the classic cprw, i.e. it lost its whole advantage, while standard wells were unaffected because one block per well makes the two conventions coincide. Sum the C and D column contractions over all of a well's block rows. That is the Galerkin column for a prolongation spreading a well's coarse unknown over all of its segment pressures by a constant, so the prolongation now spreads to match and the coarse operator stays R*S*P. The same case now takes 2494 linear iterations in 152.0 s against 2716 in 156.4 s for cprw. Default well_transfer becomes no_prolongation: the restriction is kept, the segment pressure correction is discarded and the trailing well solve corrects the wells. Prolonging one coarse value onto segment pressures as well as rates and compositions is not well defined, and it measures worse -- 2685 with the prolongation against 2494 without, while dropping the restriction too costs ~3 % (2565). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit 496fb6d)

MultisegmentWellEquations::extractCPRPressureMatrix computes one weight vector per well, averaged over every perforation of every segment, and applies it to all of the well's rows. cellavg now does the same. The previous per-block-row averaging is kept as cellblockavg. On Norne with all wells converted to multisegment the two are close -- 2547 against 2565 linear iterations -- and both beat the classic cprw's 2716, so this is not why the system solver and cprw disagree: both still diverge from cprw at the sixth linear solve rather than in the roundoff regime. The remaining difference is the coarse diagonal, which classic builds as the negated row sum without touching D while this contracts D. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit cfde17b)

StandardWellEquations::extractCPRPressureMatrix gives a pressure-controlled well a unit diagonal and skips its B and C contributions, so its coarse equation is dp = 0 rather than a contracted well equation. The system CPRW stage contracted the matrix regardless of control mode. Add preconditioner.well_identity_on_pressure_control, on by default, doing the same. Which wells are pressure controlled is decided in the outer layer and handed down as one flag per well on the WellDofLayout, so the preconditioner still sees nothing but plain data. On SPE1 this takes the system solver from 459 to 440 linear iterations against 443 for the classic cprw. It does not, however, explain why the two diverge: the per-iteration residual norms agree to six digits for the first eight and then differ by ~5e-4, at a solve where no well is pressure controlled, and that divergence is unchanged by this option, by fixing the AMG setup with --cpr-reuse-setup=3, and by using ILU0 rather than UMFPack for the well solve. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit 35dfe3b)

Classic CPRW uses two conventions: StandardWellEquations contracts D for the coarse diagonal, while MultisegmentWellEquations sets it to minus the sum of the well row's reservoir entries and never reads D. The system stage always contracted D, which made the multisegment coarse matrix differ from the very first preconditioner application. Add preconditioner.well_coarse_diagonal with auto (contract D for single-block wells, row sum for multisegment ones, i.e. classic), contract_d and row_sum. Also dump the coarse matrix and right-hand side when verbosity exceeds 10, mirroring the classic path, so the two coarse systems can be diffed directly. On SPE1 with every well converted to multisegment this brings the system solver to 436 linear iterations, exactly the classic cprw count, with the per-iteration residual norms agreeing for ten lines and first differing at line 11 -- the same point at which the well-free cpr pair diverges. The multisegment formulation difference is therefore gone and only roundoff remains. Default stays contract_d pending a decision on which is wanted generally. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit 15a984a)

ISTLSolverSystem::solve() overrides ISTLSolver::solve() and dropped the verbosity > 10 writeSystem block along with it, so --linear-solver=system_cpr silently wrote no reports/ dump. Restore it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit 397825c)

A settings sweep on full Norne, one factor at a time off the shipped defaults, serial, at the deck's own outer solver settings (maxiter 20, tol 5e-3): standard wells one segment per connection cprw 2262 2716 system_cpr 4303 4839 system_cprw, shipped defaults 2247 2569 well_transfer = classic 2248 2558 well_transfer = full 2246 2576 well_coarse_diagonal = row_sum 2242 2580 well_weight_type = cellblockavg 2247 2604 well_weight_type = quasiimpes 2245 5728 well_identity_on_pressure_control=0 2301 2574 Only well_transfer is worth changing, and only for multisegment wells: 2558 against 2569. Everything else is either flat (standard wells span 2242-2248) or worse. The margin is thin and comes from a case with an exact well solve, which nearly annihilates the well residual there is to restrict -- the other two modes stay, and are the ones to try if the well solve is ever made inexact. cellblockavg and quasiimpes are the two that must not become the default: 2604 and 5728 against 2569. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit bd08be1)

…rflow area The code fallback was "full" while setupPropertyTree ships "classic", so a JSON omitting the key silently selected a different preconditioner. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

--linear-solver-verbosity drives the Krylov output and the full-system dump; the coarse dump is a developer aid reached by setting preconditioner.verbosity in a JSON config. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

ISTLSolverSystem and SystemCprwPressureStage each had their own lookup from merged well block row to well index. Both only need the offsets, so the lookup now lives on WellDofLayout. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

buildWellDofLayout() read wellLayout_.identityOnPressureControl before prepareSystemSolver() had refreshed it from the property tree, so a pressure-controlled well got the wrong (contracted rather than trivial) coarse equation on the first linear solve of every run. Read the option first. validateSystemCPRTree checked coarsesolver/well_solver settings but not well_transfer/well_coarse_diagonal, so a typo in either only surfaced as an exception on the first linear solve, deep into a run, instead of at setup time. well_weight_type had no validation at all and silently fell through to the quasi-impes branch. Validate all three up front. buildCoarsePattern could divide by zero (NaN cast to std::size_t) on an MPI rank that locally owns no reservoir cells; guard it. dumpCoarseMatrix/dumpCoarseRhs used a rank-less filename, so MPI ranks sharing a working directory clobbered each other's debug dumps. Add the rank to the name. assembleCoarseMatrix's D-loop looked the well index back up per nonzero via wellOfBlock() even though it's already the loop variable j, and wrote through a coarse-matrix entry buildCoarsePattern never allocates for cross-well columns. Accumulate straight into the diagonal instead, which also removes that latent pattern/assembly mismatch. Also fold in outstanding style nits: reorder/alphabetize SystemCprwPressureStage.hpp's includes, rewrap lines over the actual 100-column clang-format limit (CONTRIBUTING.md's "120" is stale), and add the trailing newline three of the new JSON test fixtures were missing. Verified on full Norne: system_cprw's standard-well iteration count is unchanged at 2248, exactly matching the pre-fix baseline, and both unit test suites still pass

A pressure-controlled well gets an identity coarse row, meaning dp_j = 0, but its coarse right-hand side still received the restricted well residual. The coarse solve then returned that residual as dp_j, and full transfer prolonged it onto every segment pressure of the well. Skip these wells in the restriction, and test it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

- cellavg gives every block of a well the same weights, so compute the average once per well instead of once per block; a multisegment well with a segment per connection was quadratic. - The pressure-control flags are indexed by position in this rank's well container, the loop addBCDMatrix also uses. Say so, and throw if the count disagrees with the extracted well matrices. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

hnil · 2026-09-21T11:42:19Z

jenkins build this please

ElyesAhmed · 2026-09-22T05:39:44Z

I approved this after numerous rounds of review and testing! @atgeirr It would be nice to do one more pass!

hnil force-pushed the feature/system-cprw branch from cfde17b to b55e890 Compare August 5, 2026 12:45

hnil added the manual:new-feature This is a new feature and should be described in the manual label Aug 5, 2026

hnil mentioned this pull request Aug 6, 2026

A general system preconditioner: the CPRW pressure stage as a two-level method #7289

Draft

ElyesAhmed reviewed Aug 10, 2026

View reviewed changes

Comment thread opm/simulators/linalg/system/SystemPreconditioner.hpp Outdated

ElyesAhmed mentioned this pull request Aug 10, 2026

Let the classic CPRW contract D for multisegment wells #7282

Open

ElyesAhmed reviewed Aug 10, 2026

View reviewed changes

Comment thread opm/simulators/linalg/system/ISTLSolverSystem.hpp

Comment thread opm/simulators/linalg/system/ISTLSolverSystem.hpp Outdated

Comment thread opm/simulators/linalg/system/SystemCprwPressureStage.hpp

hnil force-pushed the feature/system-cprw branch 3 times, most recently from 6808a34 to 7adc77a Compare August 11, 2026 12:06

hnil mentioned this pull request Aug 12, 2026

System CPRW: THP well treatment and trivial-row scaling #7316

Draft

hnil requested a review from ElyesAhmed September 15, 2026 09:03

hnil force-pushed the feature/system-cprw branch from 7adc77a to 22fb7d2 Compare September 17, 2026 14:15

hnil force-pushed the feature/system-cprw branch 2 times, most recently from a1644b0 to 45b5cdf Compare September 18, 2026 11:21

ElyesAhmed marked this pull request as ready for review September 18, 2026 11:31

ElyesAhmed requested changes Sep 18, 2026

View reviewed changes

hnil and others added 11 commits September 21, 2026 13:32

Address review: match the shipped well_transfer default, note the ove…

9aceddd

…rflow area The code fallback was "full" while setupPropertyTree ships "classic", so a JSON omitting the key silently selected a different preconditioner. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

hnil and others added 6 commits September 21, 2026 13:32

Say that the coarse dump is JSON-only, on purpose

5323257

--linear-solver-verbosity drives the Krylov output and the full-system dump; the coarse dump is a developer aid reached by setting preconditioner.verbosity in a JSON config. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Address review: one wellOfBlock, on WellDofLayout

6075cd7

ISTLSolverSystem and SystemCprwPressureStage each had their own lookup from merged well block row to well index. Both only need the offsets, so the lookup now lives on WellDofLayout. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

minor cleanup

407dfd8

hnil force-pushed the feature/system-cprw branch from 45b5cdf to 3fa5301 Compare September 21, 2026 11:41

ElyesAhmed approved these changes Sep 22, 2026

View reviewed changes

Navigation Menu

Add a CPRW pressure stage to the system solver - #7278

Add a CPRW pressure stage to the system solver#7278
hnil wants to merge 17 commits into
OPM:masterfrom
hnil:feature/system-cprw

hnil commented Aug 4, 2026 •

edited

Loading

Uh oh!

Uh oh!

ElyesAhmed left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hnil commented Sep 17, 2026

Uh oh!

hnil commented Sep 18, 2026

Uh oh!

ElyesAhmed commented Sep 18, 2026

Uh oh!

ElyesAhmed left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hnil commented Sep 21, 2026

Uh oh!

ElyesAhmed commented Sep 22, 2026

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

hnil commented Aug 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

ElyesAhmed left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hnil commented Sep 17, 2026

Uh oh!

hnil commented Sep 18, 2026

Uh oh!

ElyesAhmed commented Sep 18, 2026

Uh oh!

ElyesAhmed left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hnil commented Sep 21, 2026

Uh oh!

ElyesAhmed commented Sep 22, 2026

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

hnil commented Aug 4, 2026 •

edited

Loading