A general system preconditioner: the CPRW pressure stage as a two-level method - #7289
A general system preconditioner: the CPRW pressure stage as a two-level method#7289hnil wants to merge 22 commits into
Conversation
| // subtract from the residual through C and D. | ||
| // -------------------------------------------------------------------------- | ||
| template <class Scalar, class ResComm = Dune::Amg::SequentialInformation> | ||
| class SystemCprwStep : public SystemSweepStep<Scalar> |
There was a problem hiding this comment.
maybe need a fix as SystemCprwStep is defined but never instantiated anywhere..... The CPRW pressure stage is used as the two-level coarse solver via SystemPressureBhpTransferPolicy; the sweep only ever gets ReservoirStep/WellStep. Could we either delete it, or wire it up as a block: "cprw" sweep step if the intent was to allow the pressure stage inside the sweep..?
There was a problem hiding this comment.
Right, SystemCprwStep is only reachable as the two-level coarse solver via SystemPressureBhpTransferPolicy, never from the sweep. It is there so a config can name it as a step; if you would rather it not exist until something instantiates it, I can drop it.
| return {}; | ||
| } | ||
| const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt"; | ||
| return makeWeightsCalculator(prm.get("preconditioner.weight_type"s, "quasiimpes"s), |
There was a problem hiding this comment.
This extraction is shared by every cpr/cprw configuration right?, not just the new solver, so a regression here would be wide.
There was a problem hiding this comment.
Correct - it is on the shared cpr/cprw path, so a regression would be wide. That is why the extraction is a pure refactor with the same arithmetic; worth a reviewer eye on exactly that though.
| // describe the coarse space and are ignored by FlexibleSolver. | ||
| const auto& coarsePrm = *coarseResPrm_; | ||
| const auto wellTransfer = wellTransferFromString( | ||
| coarsePrm.get("well_transfer", std::string{"full"})); |
There was a problem hiding this comment.
as before, we have a "minor" same consistency as in #7278: the code fallback here is "full" while the generated config ships "classic", so a general_system_cprw JSON omitting well_transfer silently gets a different preconditioner than the built-in one. Consider aligning the fallback to "classic" in both PRs.
There was a problem hiding this comment.
Fixed in both places here - the fallback is classic now.
ElyesAhmed
left a comment
There was a problem hiding this comment.
The current unit test does not exercise the parallel (MPI) path — that coverage is provided only by the deck runs mentioned in the commit message. This is acceptable for a unit test in isolation.
However, if you would like to strengthen the testing, consider adding a small 2-rank test case. This would specifically validate the copyOwnerToAll placement inside apply and catch any regressions in the parallel logic.
Not necessary maybe in this PR
ElyesAhmed
left a comment
There was a problem hiding this comment.
This PR turns the fixed three‑stage preconditioner into a true two‑level method (TwoLevelMethodCpr + SystemPressureBhpTransferPolicy coarse space, with an ordered sweep as the fine smoother), while exactly reproducing the original behaviour. The reproduction is structural, not accidental — the setupSystem* helpers are shared between both layouts, and the test proves equality with a well‑reasoned tolerance exception for well_transfer=full. One thing is related to SystemCprwStep which is dead code see comment).
Maybe a rebase on #7278 once it merged is a good step, so this diff drops to its own 9 commits.
27835ba to
755bdff
Compare
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)
Each step reads the residual, corrects one block, and subtracts that correction from both halves of the residual, so a step is appendable to any other and composing them is a multiplicative sweep. SystemSweepPreconditioner runs an ordered list of them as one Dune preconditioner, carrying its defect internally -- which is what a two-level method needs, since Dune hands a smoother a single (correction, defect) pair. The residual is maintained incrementally rather than recomputed as d - S v. The two are equal in exact arithmetic but not bit for bit, and the incremental form is the one the existing three-stage SystemPreconditioner uses. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit 51d443f)
pre() resizes u_ and rhs_, which a MultiTypeBlockVector cannot do, and pre() is a virtual override so it is instantiated whether or not it is called. That kept the two-level method from being used on a block system at all. Guarded with if constexpr. apply() assigns to both vectors, which sizes them, so nothing is lost for a fine level that cannot resize. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit a432593)
SystemPressureBhpTransferPolicy presents the CPRW pressure system, of dimension
nCells + nWells, as a Dune level transfer policy; PressureSolverPolicy builds
its solver; TwoLevelMethodCpr drives them with preSteps = 0 and postSteps = 1.
Everything after the coarse correction -- the coarse solver's own well solve and
then the smoother's parts -- is composed into one SystemSweepPreconditioner as
the fine smoother, so the sequence applied is
coarse -> well -> reservoir smoother -> well
which is what the fixed three-stage SystemPreconditioner does. The parts are
named by the property tree,
coarse_solver { reservoir_solver, well_solver }
smoother { reservoir_smoother, well_solver }
and every sub-tree is optional, so parts can be left out -- which the fixed
sequence cannot express. Without add_wells there is no system-wide coarse space:
the reservoir solve is an ordinary block step in the sweep and no two-level
method is built.
updateForChangedWellStructure rebuilds everything rather than refreshing in
place. That is not cosmetic: refreshing a CPR reservoir solve keeps its existing
hierarchy while rebuilding re-aggregates it, and the two are different
preconditioners. Rebuilding is the conservative choice for now.
Registered as preconditioner type general_system_cpr on the same three
operator/communication combinations as system_cpr.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
(cherry picked from commit 4ca0973)
setupGeneralSystemCPR writes the same sub-solvers setupSystemCPR does, just under coarse_solver/ and smoother/ rather than side by side, so the two configurations are the same algorithm with the same settings. The sub-solvers are now written by shared helpers, so they cannot drift apart. The JSON validation for the general layout only checks the composition, since which parts are present is the point of it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit a5792b9)
test_GeneralSystemPreconditioner applies both to the same right-hand side on a
system with a standard and a multisegment well. For well_transfer classic and
no_prolongation, and for the reservoir-only pressure stage, it compares with
BOOST_CHECK_EQUAL rather than a tolerance: TwoLevelMethodCpr reduces the defect
with applyscaleadd(-1, lhs, rhs), and SystemMatrix::usmv accumulates A,C into
the reservoir half and B,D into the well half -- the same operations in the same
order as the fixed version's four mmv calls, so the two agree exactly.
well_transfer = full is the exception and gets a tolerance. There the well half
is a sum of three non-zero terms and the two associate it differently,
((0 + dw) + c1) + c2 against dw + ((0 + c1) + c2), because Dune adds the coarse
correction to the update and then adds the smoother's separately accumulated
correction. The reservoir half has two terms and stays exact; the other transfer
modes stay exact because the coarse well correction is zero there.
Dropping either well solve changes the result, which is what makes the equality
meaningful.
Measured on decks, general against fixed, linear iterations:
system_cprw system_cpr
SPE1CASE1 440 = 440 571 = 571
SPE9_CP 468 = 468 682 = 682
Norne 2248 = 2248 4303 vs 4253
Norne msw 2558 = 2558 4839 vs 4972
with all 2831 physical summary vectors identical on Norne for the cprw pairs
(TCPU is the only one that differs) and the SPE files identical byte for byte.
The system_cpr column differs on Norne alone, the only deck here that opens and
closes wells: that path has no coarse space, so its reservoir CPR solve is
rebuilt rather than refreshed. Refreshing in place reproduces 4303 and 4839.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
(cherry picked from commit 8e873e0)
A step subtracted its correction from both halves of the residual whether or not anything read them afterwards. For a reservoir step that is a full A*corr over every cell, and in the default composition -- coarse, well, reservoir smoother, well -- nothing reads the reservoir half after the smoother. The hand-written three-stage preconditioner leaves the same updates out by hand. Which halves are still read is decided once, walking the step list backwards: a step needs the reservoir half only if a later step is a reservoir step, the well half only if a later step is a well step, and the last step needs neither. Skipping an update nothing reads cannot change a result, and the equality against the fixed preconditioner still holds. On full Norne it takes the general variant's overhead from 3.8% to about 1.7%, which is within the run-to-run spread. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit e969719)
Two things were conflated. What the change was:
Values same wells, same patterns, new numbers
Pattern same dimensions, different sparsity -- a connection opened or
closed inside a well that was already there
Dimension a well or segment appeared or vanished, so D and the coarse
system change size
and what to do about it, now preconditioner.well_structure_update:
rebuild build every part again from the property tree (default)
refresh create only what the wells size -- the well solves and the coarse
system -- and refresh the rest in place
These are different preconditioners, not two routes to one: a CPR reservoir
solve keeps its hierarchy when refreshed and re-aggregates it when rebuilt. On
full Norne, general_system_cpr gives 4253 linear iterations rebuilding and 4303
refreshing, and 4303 is what the fixed three-stage system_cpr gives. With
refresh both variants reproduce it exactly -- all 2831 physical summary vectors
identical, TCPU aside -- at 143s against 142.5s for cprw and 172s against 175.5s
for cpr.
WellMatrixStructure::hasSameDimensions separates Pattern from Dimension and
ISTLSolverSystem passes the verdict down. Both still force the coarse system to
be built again, since a changed sparsity invalidates its pattern as surely as a
changed dimension does; the distinction is carried so that a cheaper Pattern
path can be added without re-deriving it.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
(cherry picked from commit ac602cc)
A JSON array parses to a node whose children all carry an empty key, so they cannot be reached by name, and get_child_items_as_vector reads scalars rather than sub trees. get_child_list and put_child_list handle the object-array case, which a configuration describing an ordered list of steps needs. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> (cherry picked from commit 87bcc1e)
The composition was fixed in four named slots, so the one thing "general" ought
to mean could not be said in JSON. The smoother is now an ordered list:
"preconditioner": {
"type": "general_system_cpr",
"weight_type": "trueimpes",
"pre_smooth": 0, "post_smooth": 1,
"well_structure_update": "rebuild",
"coarsesolver": {
"type": "cprw_pressure",
"well_weight_type": ..., "well_transfer": ...,
"well_coarse_diagonal": ..., "well_identity_on_pressure_control": ...,
"maxiter": 1, "tol": 0.1, "solver": "loopsolver",
"preconditioner": { "type": "amg", ... }
},
"finesmoother": {
"steps": [ { "block": "well", ... },
{ "block": "reservoir", ... },
{ "block": "well", ... } ]
}
}
so steps can be reordered, repeated or left out. The names follow cpr --
coarsesolver, finesmoother, weight_type, pre_smooth, post_smooth in the same
places -- and the sub-trees are ordinary solver specs, as they were.
Three things this removes. add_wells nested inside a reservoir solver's
preconditioner was deciding whether a two-level method got built at all; that is
now coarsesolver.type, and leaving coarsesolver out means no coarse space. The
coarse pressure solver was buried at
coarse_solver.reservoir_solver.preconditioner.coarsesolver, a coarsesolver
inside a reservoir_solver that was itself the coarse solver; it is now the
coarsesolver node directly. And the keys describing the coarse space sat at the
top of the preconditioner while being read by ISTLSolverSystem rather than the
preconditioner; they now sit beside the coarse solver they belong to, and both
readers look in one place.
SPE1CASE1 440/571 and SPE9_CP 468/682, matching system_cprw and system_cpr as
before.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
(cherry picked from commit e79a9e9)
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The general layout names the weight type at the top of the preconditioner rather than in a nested CPR sub-tree. Rather than give ISTLSolver a second entry point for that, the derived solver hands the base class the two keys its existing public getWeightsCalculator() reads, so the shared cpr/cprw path is byte-identical to master. Also pass the solver verbosity down: the pressure stage dumps its coarse matrix above 10, but the sub-tree carried a hardcoded 0. 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. Reverts the config-side change from the previous commit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
6252634 to
ef7584a
Compare
Builds on #7278 — its ten commits are the first ten here, since GitHub cannot
base a pull request on a fork branch. Review from
Add composable sweep steps for the system preconditioneronwards.SystemPreconditionerruns a fixed three-stage sequence. Written out, thatsequence is
which is
coarse_solver { reservoir_solver, well_solver }followed bysmoother { reservoir_smoother, well_solver }. This PR makes that structure theconfiguration, and makes the pressure stage a real two-level method rather than a
hand-rolled first stage.
SystemPressureBhpTransferPolicypresents the CPRW pressure system (dimensionnCells + nWells) as a Dune level transfer policy,
PressureSolverPolicybuildsits solver, and
TwoLevelMethodCprdrives them withpreSteps = 0,postSteps = 1. Dune hands a smoother a single (correction, defect) pair, soeverything after the coarse correction is composed into one
SystemSweepPreconditionerthat carries its own defect. Withoutadd_wellsthere is no system-wide coarse space, so the reservoir solve is an ordinary block
step and no two-level method is built.
Selected with
--linear-solver=general_system_cpr/general_system_cprw. TheJSON follows cpr's shape --
coarsesolver,finesmoother,weight_type,pre_smooth,post_smoothin the same places -- but the smoother is an orderedlist, since a coupled system has more than one block to visit:
so steps can be reordered, repeated or left out.
coarsesolveris itself thesolver spec for the coarse pressure system and the keys beside it describe the
coarse space; leaving it out means no coarse space. That replaces an
add_wellsburied in a nested preconditioner deciding whether a two-level method got built
at all, and a coarse solver reached at
coarse_solver.reservoir_solver.preconditioner.coarsesolver.Arrays of sub-trees needed
PropertyTree::get_child_list/put_child_list,which did not exist --
get_child_items_as_vectorreads scalars.With equivalent settings it reproduces the existing preconditioner exactly.
That works because
TwoLevelMethodCprkeeps the same defect bookkeeping thesweep does --
postsmoothreduces the defect withapplyscaleadd(-1, lhs, rhs),and
SystemMatrix::usmvaccumulates A,C into the reservoir half and B,D into thewell half, the same operations in the same order as the fixed version's four
mmvcalls.system_cprwgeneral_system_cprwsystem_cprgeneral_system_cprOn Norne all 2831 physical summary vectors are identical for both
system_cprwpairs; TCPU is the only one that differs. The two
general_system_cprfiguresthat do not match are the
updateForChangedWellStructuretrade described below --Norne is the only one of these decks that opens and closes wells.
One exception, which the test states explicitly: with
well_transfer = fullthewell half of the result is a sum of three non-zero terms, and the two associate it
differently --
((0 + dw) + c1) + c2againstdw + ((0 + c1) + c2), because Duneadds the coarse correction to the update and then adds the smoother's separately
accumulated correction. They agree to one rounding step and no closer. The
reservoir half has two terms and stays exact, and
classic(the default) andno_prolongationstay exact because the coarse well correction is zero there.Removing this would mean rewriting Dune's
postsmooth, which also serves theclassic CPR path.
Two incidental changes:
TwoLevelMethodCpr::preresizes its work vectors, andMultiTypeBlockVectorhas no
resize. Sincepreis a virtual override it is instantiated whether ornot it is called, so the resize is now behind
if constexpr (requires ...).applyassigns to both vectors, which sizes them.preconditioner.well_structure_updatechooses what a changed well structurecosts:
rebuild(default) builds every part again,refreshcreates only whatthe wells size and refreshes the rest in place. Not cosmetic -- a CPR reservoir
solve keeps its hierarchy when refreshed and re-aggregates it when rebuilt, and
on Norne that moves
general_system_cprbetween 4253 and 4303 iterations.refreshis what the fixed version does, and with it both variants reproduce itexactly. The kind of change is classified too (
Values/Pattern/Dimension), thoughPatternandDimensionare still handled alike: achanged sparsity invalidates the coarse pattern as surely as a changed dimension.
Cost. Full Norne, two repeats on one binary, mean wall time:
rebuildrefreshsystem_cprwsystem_cprRun-to-run spread is 3-5%, so nothing here separates them. The first version was
a consistent 3.8% slower because each step updated both halves of the residual
whether or not anything read them again -- a full
A*corrover every cell in thesmoother, which the hand-written version omits by hand. Which halves are still
read is now decided once at build time.
🤖 Generated with Claude Code