System CPRW: THP well treatment and trivial-row scaling - #7316
Draft
hnil wants to merge 28 commits into
Draft
Conversation
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>
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>
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>
First half of the general two-level infrastructure, alongside the existing 3-stage preconditioner rather than replacing it. SystemPreconditionerParts provides the pieces a coupled-system preconditioner is built from: SystemWellSolve corrects the well unknowns from the defect d[_1] - B v[_0] - D v[_1], SystemReservoirSolve does the same for the reservoir block, and SystemReservoirWellSweep pairs them into one multiplicative sweep -- a reservoir smoother followed by a well solve. Each part recomputes its own defect, so composing them is block Gauss-Seidel rather than block Jacobi. SystemPressureBhpTransferPolicy presents the CPRW pressure system as a Dune LevelTransferPolicyCpr so TwoLevelMethodCpr can drive it, delegating the assembly and the transfers to SystemCprwPressureStage rather than repeating them. clone() forces the policy to be copyable while the stage owns unique pointers, so the stage is held by shared_ptr and clone() is a shallow copy, as it is for PressureBhpTransferPolicy. The stage grows caller-supplied-output forms of the two transfers so the policy needs no extra copies, and both now share one implementation. 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>
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>
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>
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>
verbosity above 10 now writes the coarse pressure matrix from the classic PressureBhpTransferPolicy as well as from the system stage, so the two can be diffed entry by entry instead of inferring differences from residual norms several steps downstream. On SPE1 with standard wells this localises the disagreement: the two coarse matrices are bit identical for the first four assemblies and then differ by about 10 % in a reservoir-row, well-column entry -- the C contraction for one well -- after which the difference spreads into the reservoir entries. That is not roundoff. The leading suspect is when the two ask whether a well is pressure controlled: classic asks while assembling inside the preconditioner update, this asks once per prepareSystemSolver in buildWellDofLayout, so a control switch between those points makes one zero the column and the other not. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Adds a weights dump next to the coarse matrix dump, under the same verbosity > 10 switch, on both the classic transfer policy and the system stage. This locates the standard-well disagreement between cprw and system_cprw. On SPE1 the weight vectors are bit identical for the first four assemblies and then differ by a factor of nearly two -- 0.00407767 against 0.0020695 at the last component of the last cell. The coarse rows for that cell then differ by a uniform relative 9.776e-02, which is what a differing weight vector does to a row, and the residuals follow. So the disagreement is not in the well coupling at all: not the C contraction, the coarse diagonal, the control mode or the well solver, all of which were checked and cleared. Both sides build the same trueimpes calculator over the same element context, so the remaining explanation is that one of them evaluates it against a different simulator state, i.e. a weight vector that is not refreshed when it should be. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
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>
Diagnostic, verbosity > 10, next to the coarse-matrix and weights dumps: primvars_N.txt holds every cell's primary variables and the gas/water/pressure meanings. Written from initPrepare(), which both the classic and the system solver path go through, so the two can be compared cell by cell. This is what settled the cprw vs system_cprw disagreement at assembly 4: the states agree to 1e-15 and 17 cells have simply picked the other side of the Sg/Rs switch, which changes the trueimpes weights without changing the state. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
preconditioner.well_coarse_diagonal = contract_d takes a well's coarse diagonal from lambda' D(:,p) instead of minus the row sum of its contracted reservoir entries. Default "auto" keeps today's behaviour. Standard wells have always contracted D (StandardWellEquations reads duneD_[0][0]); only the multisegment path used the row sum, which never reads D and therefore throws away all segment-to-segment coupling. With one segment per connection that is most of the well. Contracting D is also the Galerkin diagonal for the prolongation the coarse column already assumes -- one coarse value spread over all of the well's segment pressures. The flag is threaded next to use_well_weights, from PressureBhpTransferPolicy down to MultisegmentWellEquations::extractCPRPressureMatrix. The contraction itself is mswellhelpers::contractCprWellDiagonal, unit tested in test_MswCprWellDiagonal; it returns 1 on exact cancellation, since a zero would make the coarse pressure system singular. Norne with --convert-to-multisegment-well=per-connection: every coarse matrix differs from the default, i.e. the flag reaches the assembly. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
The classic CPRW gives every pressure-controlled well a trivial coarse row, and isPressureControlled counts THP as pressure control. The trivial row assumes constant bhp; the contracted row assumes constant rates; under THP control neither holds, since bhp and rates co-vary along the VFP curve. The new key (default true, matching classic) lets a THP well keep its contracted equation, reserving the trivial row for actual BHP control, where it is exact. On model5/5_NETWORK_MODEL5_STDW (network, non-monotonic VFP, wells oscillating ORAT<->THP) the contracted THP row takes 232 linear iterations against 395 for the trivial row on an identical Newton path; the MSW variant gives 306/195 against 406/204. Also reads the coarse-space options before building the well DOF layout, so the first assembly of a run honours them; previously the first Newton iteration used the defaults. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Extends the cellavg well weights with a control-row weight mu = -(w'D_cq)/D_kq for THP-controlled standard wells, chosen to cancel the WQTotal column of the contracted equation. The coarse diagonal becomes the Schur value w'D_cb - (w'D_cq) D_kb/D_kq, i.e. the total-rate correction is eliminated through the control equation instead of being frozen. Since B's control row is zero the reservoir couplings are unchanged, and rate control has D_kb = 0, so only THP wells are affected. A flat VFP (D_kq ~ 0) degenerates to BHP control and keeps the trivial row; a sign-flipping or collapsing corrected diagonal falls back to the plain contracted row; multisegment wells are excluded. All of it lives in the outer layer. Verified against dumps on model5/5_NETWORK_MODEL5_STDW: the only coarse-matrix entry that changes is the THP well's diagonal, and it changes by exactly mu. At a tight linear tolerance all variants share one Newton path and measure within 3%: contracted 1346, trivial 1368, cellavg_vfp 1390, classic cprw 1431. The Schur diagonal does not pay on this deck; the option exists to test decks where THP wells are strongly rate-responsive. Also dumps the well part of the CPR weights next to the reservoir part at verbosity > 10. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Replaces well_thp_is_pressure_control and the cellavg_vfp weight type with one key saying how a THP-controlled well enters the coarse system: bhp (trivial row, classic behaviour, default), rate (keep the contracted equation), or vfp (contracted equation plus the control-row elimination). BHP wells always get the trivial row and rate-controlled wells the contracted equation regardless, since those cases are not in question. Verified to reproduce the previous measurements exactly on model5/5_NETWORK_MODEL5_STDW at tight tolerance: 1368 / 1346 / 1390. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A pressure-controlled well's trivial coarse row is decoupled and has a zero right-hand side, so its diagonal cannot affect the coarse solution - but it leaks into the AMG anyway. The coarse pattern keeps the well couplings as stored zeros, the aggregation follows the pattern, and the Galerkin product mixes the diagonal into whatever aggregate the well node joins. The previous value of 1.0 sits ~1e6 above the reservoir coarse diagonals and poisons the coarse correction near every pressure-controlled well. Measured on model5/5_NETWORK_MODEL5_STDW (network, wells oscillating between ORAT and THP): a scale sweep of the trivial diagonal is monotone and saturates in both directions (1e3: 1365, 1.0: 1368, 1e-3: 1365, reservoir scale: 1319, 1e-12: 1224 linear iterations at tight tolerance; at production tolerance the default drops 395 -> 184 and takes a better Newton path). Norne is unaffected (2248 -> 2242) - its history-controlled wells are essentially never pressure-controlled. The fix sets the trivial diagonal to 1e-6 times the smallest reservoir coarse diagonal, i.e. well into the saturated harmless regime. With the scale artifact removed, the THP-treatment comparison inverts: the trivial row measures best (1223), ahead of the contracted row (1346) and the control-row elimination (1317). Also zero the restricted right-hand side of pressure-controlled wells under the full and no_prolongation transfers; restricting a residual against the trivial equation would produce an arbitrary bhp correction, which the small diagonal would amplify. The same 1.0 sits in the classic CPRW (StandardWellEquations:: extractCPRPressureMatrix, whose comment already suspected the scaling) over the same stored-zero pattern; fixing that is a separate, upstream change. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Same defect and same fix as the system CPRW stage: the trivial equation of a pressure-controlled well is decoupled with a zero right-hand side, so its diagonal cannot affect the solution, but the coarse pattern keeps the couplings as stored zeros, the AMG aggregation follows the pattern, and the old value of 1.0 - far above the reservoir rows' scale - poisons whatever aggregate the well node joins. The comment at the standard-well diagonal already suspected the scaling. The cell-average weights are now accumulated for pressure-controlled wells as well (couplings still not written), and the trivial diagonal becomes 1e-6 times the contracted diagonal's magnitude. The multisegment path additionally gains a guard for a well with no local perforations, which previously divided by zero past an assert. Measured on model5/5_NETWORK_MODEL5_STDW with --linear-solver=cprw: 1431 -> 1289 linear iterations at linear reduction 1e-8 on an identical Newton path; 308 -> 243 at the production tolerance; the MSW variant 279 -> 257. Norne is unchanged (2262 -> 2244), as its history-controlled wells are essentially never pressure-controlled. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
hnil
added
the
manual:enhancement
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Follow-up to #7278 — merge after it (and after the
general_system_cprcommits it builds on; only the last 5 commits here are new).Three things:
well_thp_treatment=bhp|rate|vfp(system CPRW) — how a THP-controlled well enters the coarse system. The trivial row assumes constant bhp, the contracted row constant rates; under THP neither holds, andvfpeliminates the rate correction through the control equation, so the coarse diagonal becomes the Schur valuew'D_cb - (w'D_cq) D_kb/D_kqcarrying the VFP slope. Defaultbhp= classic behaviour. BHP wells always get the trivial row (exact); rate-controlled wells the contracted equation.Trivial-row diagonal scaled to the problem (system CPRW stage). The old value 1.0 sits ~1e6 above the reservoir coarse diagonals; the pattern keeps the well couplings as stored zeros, AMG aggregation follows the pattern, and the Galerkin product mixes the off-scale diagonal into the aggregate — degrading the coarse correction near every pressure-controlled well. A scale sweep is monotone and saturates both ways.
The same fix in the classic CPRW (
StandardWellEquations/MultisegmentWellEquations, whosediagElem = 1.0comment already suspected the scaling): the trivial diagonal becomes1e-6 * |contracted diagonal|.Measured on
model5/5_NETWORK_MODEL5_STDW(network, non-monotonic VFP, wells oscillating ORAT/THP), linear iterations. "Tight" = linear reduction 1e-8, which puts every solver on the identical Newton path so the counts compare preconditioners, not trajectories; "default" = linear reduction 0.005. Nonlinear tolerances untouched.cprwbeforecprwaftersystem_cprwafterMSW variant: cprw 279 → 257, system 258. SPE1 443 → 435, SPE9 439 → 430 (system vs classic-before), Norne serial cprw 2262 → 2244, system 2242 (unchanged — its wells are essentially never pressure-controlled). THP treatments at tight reduction:
bhp1223,vfp1317,rate1346;vfpcan destabilise the Newton path on control-oscillating decks at the default reduction, so it stays an expert option.🤖 Generated with Claude Code