Skip to content

Add optional scaling of the well bhp and rate primary variables - #7286

Open
hnil wants to merge 3 commits into
OPM:masterfrom
hnil:pr/well-primary-variable-scaling
Open

Add optional scaling of the well bhp and rate primary variables#7286
hnil wants to merge 3 commits into
OPM:masterfrom
hnil:pr/well-primary-variable-scaling

Conversation

@hnil

@hnil hnil commented Aug 6, 2026

Copy link
Copy Markdown
Member

The bhp column of the well D block sits ~7 decades below the rate column because bhp is in Pascals: measured cond(D) is 1e9 (SPE1) to 1e10 (thermal). This adds --well-bhp-scaling and --well-rate-scaling, applied to the AD derivative only, so stored values stay physical and B, C, D and the Schur complement remain consistent; C D^-1 B is invariant and iteration counts are unchanged.

Defaults are 1.0 and bit-identical off. Recommended setting: --well-bhp-scaling=8388608 (2^23, ~84 bar), which brings cond(D) to ~4e3 and moves the smallest |det| from 3e-11 to 2e-4 — seven decades further from the absolute 1e-40 singularity threshold. Rate scaling measures as already O(1), hence its default. Main gain is robustness for float builds and threshold tests; open question is whether the default should become 2^23 after wider testing.

Shows the effect (scaling active, results and iteration counts unchanged):

flow SPE1CASE1.DATA --linear-solver=cprw --well-bhp-scaling=8388608 --well-rate-scaling=4

Covers standard and multisegment wells. Verified at defaults against SPE1CASE1 (437), BASE2_MSW_HFA (30), MSW-2D-HZ (356), SPE1CASE2_MSW_THERMAL (73) and SPE1CASE2_THERMAL (72).

Depends on #7314, which fixes a hard-coded Jacobian entry that a non-unit --well-rate-scaling would otherwise expose.

On the representation: an alternative is to store the scaled variable x/s and multiply out in the accessors, which keeps the stored value and the linear system in one space. This PR keeps the stored value physical instead, so the conversion is confined to the Newton increment (one helper per class) rather than every WellState exchange and getter, and restart/log content stays independent of the parameter. Both are documented next to value_.

🤖 Generated with Claude Code

@hnil

hnil commented Aug 6, 2026

Copy link
Copy Markdown
Member Author

Extended to multisegment wells (same derivative-only design; segment pressure and total rate). One consumer needed a matching fix: the classic CPRW row-sum well diagonal estimates the physical bhp derivative and must carry the scale, otherwise linear iterations double on MSW decks. With it, scaling is iteration-neutral also for MSW (BASE2_MSW_HFA 30/30, MSW-2D-HZ 356/356, SPE1CASE2_MSW_THERMAL 73/73 with and without 2^23).

@hnil hnil added the manual:enhancement This is an enhancement/improvent that needs to be documented in the manual label Aug 7, 2026
@hnil
hnil marked this pull request as draft August 7, 2026 15:35
@hnil
hnil marked this pull request as ready for review August 10, 2026 09:34
@hnil
hnil requested a review from GitPaean August 10, 2026 09:34
@bska

bska commented Aug 10, 2026

Copy link
Copy Markdown
Member

jenkins build this please

@hnil
hnil force-pushed the pr/well-primary-variable-scaling branch from bb207e2 to b21c843 Compare August 11, 2026 08:23
@GitPaean

Copy link
Copy Markdown
Member

can you share some test findings for cases, we probably want to incorporate this in some regression tests by adding running arguments.

hnil and others added 2 commits August 11, 2026 14:06
The bhp column of the well D block sits 6-7 decades below the others: the
conservation-equation derivatives w.r.t. bhp are ~1e-7 because bhp is in
Pascals, while the control equation contributes d(bhp)/d(bhp) = 1. Measured
cond(D) medians are 9.7e8 (SPE1CASE1) and 1.3e10 (SPE1CASE2_THERMAL), fully
removable by diagonal scaling - i.e. units, not physics.

Scale only the derivative when the well Evaluations are created:

    createVariable(totalNumEq, value_[eqIdx]/s, numEq+eqIdx) * s

The stored value_ stays physical, so update(), copyToWellState(), the
absolute bhp lower limit, the convergence checks and getPrimaryVars are all
untouched; only the Newton increments convert back (they arrive in scaled
units). B, C, D and resWell_ stay mutually consistent because they are all
assembled from the same scaled Evaluations, and C D^-1 B is exactly
invariant under the column scaling - confirmed: iteration counts are
unchanged (437 -> 437, 72 -> 73).

--well-bhp-scaling=8388608 (2^23, ~84 bar) equilibrates: cond(D) median
drops to 3.8e3 / 3.7e3 and min |det| moves from 2.9e-11 to 2.4e-4, seven
decades further from matrixblock.hh's absolute 1e-40 branch. The gain is
robustness (float builds, absolute thresholds), not iterations.
--well-rate-scaling exists for symmetry; the rate column measures as
already O(1). Defaults are 1.0 and reproduce unscaled results bit for bit.

A per-well scale derived from the current bhp was tried and rejected:
better median conditioning but ~8x worse in the tail, because stopped and
zero-rate wells are not scaled by their bhp magnitude - and the tail is
what the scaling exists for.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Same design as for standard wells: the stored values stay physical and only
the Evaluation derivative carries the scale, so B, C, D and the Schur
complement remain mutually consistent by construction. The Newton increments
for segment pressure and total rate arrive in scaled units and are converted
in updateNewton(), where the physical limits (max_pressure_change, the bhp
lower limit) then apply unchanged. Fractions and the temperature variable
keep scale 1.

One consumer needed a matching fix: the classic CPRW coarse well row builds
its column entries from C's segment-pressure column, which carries the
scale, while the default row-sum diagonal estimates the physical bhp
derivative from unscaled B entries. Without the same factor the coarse well
column and diagonal disagree by the scale, which doubled linear iterations
on BASE2_MSW_HFA (30 -> 59) and SPE1CASE2_MSW_THERMAL (73 -> 139). With it,
scaling is iteration-neutral (30/356/73 with and without 2^23), as the
C D^-1 B invariance requires. The contracted-diagonal convention needs no
fix since it inherits the scale from D itself.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@hnil
hnil force-pushed the pr/well-primary-variable-scaling branch from b21c843 to 9e1abc3 Compare August 11, 2026 12:06
The scaling puts value_ and the linear system's unknown in different spaces,
which is easy to get wrong - the trivial-equation bug fixed in the previous
commit is exactly that mistake. Make the boundary explicit instead of
implicit: physicalIncrement() replaces the bare varScale() multiplies, so
the two places solver-space quantities enter each class are named, and the
contract is stated once next to value_ (including the alternative
representation and why it was not chosen).

Extend the unit test to pin the contract rather than the representation:
eval(i).value() == value(i) and a zero off-diagonal derivative for every
slot, value()/setValue() as exact inverses (the getPrimaryVars round trip
NLDD relies on), and a Newton step with the absolute bhp floor ACTIVE - the
existing case uses a zero increment and never reaches that branch, which is
where a missing conversion on a physical limit would show up.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@hnil
hnil force-pushed the pr/well-primary-variable-scaling branch from 4bc4d7f to 027f004 Compare August 12, 2026 08:47
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

manual:enhancement This is an enhancement/improvent that needs to be documented in the manual

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants