Add optional scaling of the well bhp and rate primary variables - #7286
Open
hnil wants to merge 3 commits into
Open
Add optional scaling of the well bhp and rate primary variables#7286hnil wants to merge 3 commits into
hnil wants to merge 3 commits into
Conversation
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
added
the
manual:enhancement
hnil
marked this pull request as draft
August 7, 2026 15:35
hnil
marked this pull request as ready for review
August 10, 2026 09:34
Member
|
jenkins build this please |
hnil
force-pushed
the
pr/well-primary-variable-scaling
branch
from
August 11, 2026 08:23
bb207e2 to
b21c843
Compare
Member
|
can you share some test findings for cases, we probably want to incorporate this in some regression tests by adding running arguments. |
hnil
force-pushed
the
pr/well-primary-variable-scaling
branch
2 times, most recently
from
August 12, 2026 08:47
4bc4d7f to
027f004
Compare
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>
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
force-pushed
the
pr/well-primary-variable-scaling
branch
from
August 19, 2026 07:42
027f004 to
a2515fd
Compare
Member
Author
|
The numbers behind this, since they are the argument:
The bhp column sits ~7 decades below the rate column simply because bhp is in Pascals. 2^23 is ~84 bar, which is why it lands near unity. The last row is the one that matters for robustness: 2e-4 is seven decades further from the 1e-40 threshold the singularity check uses, so float builds and threshold tests stop being marginal. Rate scaling measures as already O(1), hence its default of 1.0. Scaling is applied to the AD derivative only, so stored values stay physical, B/C/D and the Schur complement stay consistent, and C D^-1 B is invariant - iteration counts are unchanged at defaults on SPE1CASE1 (437), BASE2_MSW_HFA (30), MSW-2D-HZ (356), SPE1CASE2_MSW_THERMAL (73) and SPE1CASE2_THERMAL (72). The open question is whether 2^23 should become the default once it has had wider testing. |
Member
Author
|
jenkins build this please |
Contributor
There was a problem hiding this comment.
Pull request overview
Adds configurable well BHP/rate primary-variable scaling to improve Jacobian conditioning while retaining physical stored values.
Changes:
- Adds BHP and rate scaling parameters.
- Applies scaling to standard and multisegment well derivatives and Newton updates.
- Updates MSW CPR extraction and adds standard-well tests.
Reviewed changes
Copilot reviewed 8 out of 8 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
tests/test_wellmodel.cpp |
Tests standard-well scaling behavior. |
opm/simulators/wells/StandardWellPrimaryVariables.hpp |
Defines standard-well scaling contract. |
opm/simulators/wells/StandardWellPrimaryVariables.cpp |
Scales derivatives and Newton increments. |
opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp |
Defines MSW scaling contract. |
opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp |
Applies scaling to segment variables. |
opm/simulators/wells/MultisegmentWellEquations.cpp |
Adjusts MSW CPR diagonal scaling. |
opm/simulators/flow/BlackoilModelParameters.hpp |
Declares scaling parameters and defaults. |
opm/simulators/flow/BlackoilModelParameters.cpp |
Registers command-line parameters. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Comment on lines
134
to
+137
| evaluation_[eqIdx] = | ||
| EvalWell::createVariable(totalNumEq, | ||
| value_[eqIdx], | ||
| Indices::numEq + eqIdx); | ||
|
|
||
| value_[eqIdx] / s, | ||
| Indices::numEq + eqIdx) * s; |
Comment on lines
+136
to
+139
| Parameters::Register<Parameters::WellBhpScaling<Scalar>> | ||
| ("Scaling of the well bhp primary variable (1 disables it, 8388608 equilibrates)"); | ||
| Parameters::Register<Parameters::WellRateScaling<Scalar>> | ||
| ("Scaling of the well total-rate primary variable"); |
Comment on lines
+188
to
+189
| //! - eval(i) has a physical value and derivative varScale(i), so B, C and D | ||
| //! carry the scaled column while the residual rows stay physical. |
Comment on lines
+227
to
+230
| // 2^16 rather than the recommended 2^23: SetDefault round-trips the value | ||
| // through text at 6 significant digits, so it must be exactly representable | ||
| // there (8388608 would arrive as 8388610). Command-line parsing is exact. | ||
| Opm::Parameters::SetDefault<Opm::Parameters::WellBhpScaling<double>>(65536.0); // 2^16 |
Comment on lines
+81
to
+84
| // The value stays physical; the derivative d(x)/d(x/s) = s makes | ||
| // this the Jacobian column of the scaled variable. Newton | ||
| // increments arrive scaled and are converted in updateNewton(). | ||
| evaluation_[seg][eq_idx].setDerivative(eq_idx + Indices::numEq, varScale(eq_idx)); |
Comment on lines
+448
to
+450
| static const Scalar bhp_scale = | ||
| Parameters::Get<Parameters::WellBhpScaling<Scalar>>(); | ||
| diag_ell *= bhp_scale; |
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
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.
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-scalingand--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):
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-scalingwould 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