From db43a611dc947d2a3677cecce44c7ee0660f2ef9 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 14:17:14 +0200 Subject: [PATCH 1/3] Add scaling of the well bhp and rate primary variables 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 --- .../flow/BlackoilModelParameters.cpp | 4 ++ .../flow/BlackoilModelParameters.hpp | 13 ++++ .../wells/StandardWellPrimaryVariables.cpp | 42 +++++++++++-- .../wells/StandardWellPrimaryVariables.hpp | 4 ++ tests/test_wellmodel.cpp | 62 +++++++++++++++++++ 5 files changed, 119 insertions(+), 6 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelParameters.cpp b/opm/simulators/flow/BlackoilModelParameters.cpp index 1e05f78a23c..c27a81a87c6 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -133,6 +133,10 @@ void BlackoilModelParameters::registerParameters() { Parameters::Register> ("Maximum relative change of the bottom-hole pressure in a single iteration"); + Parameters::Register> + ("Scaling of the well bhp primary variable (1 disables it, 8388608 equilibrates)"); + Parameters::Register> + ("Scaling of the well total-rate primary variable"); Parameters::Register> ("Maximum absolute change of a well's volume fraction in a single iteration"); Parameters::Register> diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index caadd0fc100..3b83e347cc6 100644 --- a/opm/simulators/flow/BlackoilModelParameters.hpp +++ b/opm/simulators/flow/BlackoilModelParameters.hpp @@ -93,6 +93,19 @@ struct ToleranceWells { static constexpr Scalar value = 1e-4; }; template struct ToleranceWellControl { static constexpr Scalar value = 1e-7; }; +//! \brief Scaling of the well bhp primary variable. +//! \details The bhp column of the well D block sits ~1e-7 below the rate column +//! because bhp is in Pascals; 8388608 (2^23, ~84 bar) equilibrates it and +//! removes 5-7 orders of magnitude of condition number. Powers of two keep +//! x/s then *s exact. 1.0 disables the scaling. +template +struct WellBhpScaling { static constexpr Scalar value = 1.0; }; + +//! \brief Scaling of the well total-rate primary variable. +//! \details Measures as already O(1), so 1.0 is the right default. +template +struct WellRateScaling { static constexpr Scalar value = 1.0; }; + struct MaxWelleqIter { static constexpr int value = 30; }; template diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index ed1f89853ba..bc31c5e4f2b 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp @@ -25,6 +25,11 @@ #include #include +#include +#include + +#include + #include #include @@ -95,6 +100,25 @@ Scalar relaxationFactorFraction(const Scalar old_value, namespace Opm { +template +typename FluidSystem::Scalar +StandardWellPrimaryVariables::varScale(const int eqIdx) +{ + // A per-well scale derived from the current bhp was tried and rejected: it + // improves the median conditioning of D but is ~8x worse in the tail + // (stopped and zero-rate wells are not scaled by their bhp magnitude), and + // the tail is what the scaling exists for. + if (eqIdx == WQTotal) { + static const Scalar s = Parameters::Get>(); + return s; + } + if (eqIdx == Bhp) { + static const Scalar s = Parameters::Get>(); + return s; + } + return Scalar{1}; // the fractions are dimensionless +} + template void StandardWellPrimaryVariables:: setEvaluationsFromValues() @@ -102,11 +126,15 @@ setEvaluationsFromValues() // total number of equations/derivatives (well + reservoir) const int totalNumEq = numWellEq_ + Indices::numEq; for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) { + // value_ stays in physical units; only the derivative is scaled, so the + // Jacobian column is that of X = x/s while update(), copyToWellState(), + // the absolute bhp limit and the convergence checks keep working on + // physical values. Newton increments are converted in updateNewton(). + const Scalar s = varScale(eqIdx); evaluation_[eqIdx] = EvalWell::createVariable(totalNumEq, - value_[eqIdx], - Indices::numEq + eqIdx); - + value_[eqIdx] / s, + Indices::numEq + eqIdx) * s; } } @@ -295,7 +323,8 @@ updateNewton(const BVectorWell& dwells, this->processFractions(); // updating the total rates Q_t - value_[WQTotal] -= dwells[0][WQTotal]; + // The solution is in scaled units, value_ is physical, so convert the increment. + value_[WQTotal] -= varScale(WQTotal) * dwells[0][WQTotal]; // here, we make sure it is zero for wells with zero rate target(including stopped wells) if (stop_or_zero_rate_target) { @@ -310,8 +339,9 @@ updateNewton(const BVectorWell& dwells, } // updating the bottom hole pressure - const int sign1 = dwells[0][Bhp] > 0 ? 1: -1; - const Scalar dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), + const Scalar dbhp = varScale(Bhp) * dwells[0][Bhp]; + const int sign1 = dbhp > 0 ? 1: -1; + const Scalar dx1_limited = sign1 * std::min(std::abs(dbhp), std::abs(value_[Bhp]) * dBHPLimit); // some cases might have defaulted bhp constraint of 1 bar, we use a slightly smaller value as the bhp lower limit for Newton update // so that bhp constaint can be an active control when needed. diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index 8fb7a74f5ac..6364903620c 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -154,6 +154,10 @@ class StandardWellPrimaryVariables { DeferredLogger& deferred_logger) const; private: + //! \brief Scaling of well primary variable \p eqIdx, 1.0 unless configured. + //! \details Only the Jacobian column is scaled; value_ stays physical. + static Scalar varScale(const int eqIdx); + //! \brief Initialize evaluations from values. void setEvaluationsFromValues(); diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index 70ecd438920..9dcceb2fca8 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -52,6 +52,7 @@ #include #include +#include #include #include @@ -218,3 +219,64 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { BOOST_CHECK(well->numStaticWellEq== 4); } } + +BOOST_AUTO_TEST_CASE(TestPrimaryVariableScaling) { + // The scaling must appear in the derivative only: the stored value and the + // Evaluation's value stay physical. Set the parameters before the first + // varScale() call - the scales are cached statically. + // 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>(65536.0); // 2^16 + Opm::Parameters::SetDefault>(0.25); // 2^-2 + + const SetupTest setup_test; + const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep); + const Opm::BlackoilModelParameters param; + + using FluidSystem = Opm::BlackOilFluidSystem; + // Just enough initialisation for phase-usage queries to work + // (same device as test_rftcontainer.cpp). + FluidSystem::initBegin(/*numPvtRegions=*/1); + using RateConverterType = Opm::RateConverter:: + SurfaceToReservoirVoidage>; + RateConverterType rateConverter(std::vector(10, 0)); + + const auto& well_ecl = wells_ecl[0]; // PROD1 + Opm::PerforationData dummy; + std::vector> pdata(well_ecl.getConnections().size(), dummy); + for (auto c = 0*pdata.size(); c < pdata.size(); ++c) { + pdata[c].ecl_index = c; + } + Opm::ParallelWellInfo pinfo{well_ecl.name()}; + const StandardWell well(well_ecl, pinfo, setup_test.current_timestep, + param, rateConverter, 0, 3, 3, 0, pdata); + + using PV = std::decay_t; + PV pv(well); + pv.resize(well.numStaticWellEq); + + // Zero Newton update through the public interface triggers + // setEvaluationsFromValues(); the bhp lands on its lower limit. + typename PV::BVectorWell dwells(1); + dwells[0].resize(well.numStaticWellEq); + dwells[0] = 0.0; + Opm::DeferredLogger logger; + pv.updateNewton(dwells, /*stop_or_zero_rate_target=*/false, + /*dFLimit=*/0.2, /*dBHPLimit=*/0.1, logger); + + constexpr int numEq = StandardWell::Indices::numEq; + + // Values are physical, with and without scaling. + BOOST_CHECK_EQUAL(pv.eval(PV::Bhp).value(), pv.value(PV::Bhp)); + BOOST_CHECK_EQUAL(pv.eval(PV::WQTotal).value(), pv.value(PV::WQTotal)); + + // The derivatives carry the scale: d(x)/d(x/s) = s. These fail without the + // scaling support (both were hard-coded 1.0). + BOOST_CHECK_EQUAL(pv.eval(PV::Bhp).derivative(numEq + PV::Bhp), 65536.0); + BOOST_CHECK_EQUAL(pv.eval(PV::WQTotal).derivative(numEq + PV::WQTotal), 0.25); + + // The dimensionless fractions stay unscaled. + BOOST_CHECK_EQUAL(pv.eval(PV::WFrac).derivative(numEq + PV::WFrac), 1.0); + BOOST_CHECK_EQUAL(pv.eval(PV::GFrac).derivative(numEq + PV::GFrac), 1.0); +} From 9e1abc3e86a502379f04689180d7484931f71874 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 15:17:39 +0200 Subject: [PATCH 2/3] Extend well primary-variable scaling to multisegment wells 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 --- .../wells/MultisegmentWellEquations.cpp | 11 +++++++ .../MultisegmentWellPrimaryVariables.cpp | 32 ++++++++++++++++--- .../MultisegmentWellPrimaryVariables.hpp | 4 +++ 3 files changed, 43 insertions(+), 4 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index 7e639e53ab2..b3a7a16713a 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -35,6 +35,9 @@ #include #endif +#include + +#include #include #include #include @@ -438,6 +441,14 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, } } + // The well-column entries above are derivatives w.r.t. the *scaled* + // segment pressure (see MultisegmentWellPrimaryVariables); the row-sum + // diagonal estimates the physical bhp derivative and must carry the + // same factor, or the coarse well column and diagonal disagree by it. + static const Scalar bhp_scale = + Parameters::Get>(); + diag_ell *= bhp_scale; + #define EXTRA_DEBUG_MSW 0 #if EXTRA_DEBUG_MSW if (diag_ell <= 0.0) { diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 8f441a275b2..63046daa5e1 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -31,6 +31,9 @@ #include #include +#include + +#include #include #include #include @@ -52,6 +55,21 @@ resize(const int numSegments) evaluation_.resize(numSegments); } +template +typename FluidSystem::Scalar +MultisegmentWellPrimaryVariables::varScale(const int eqIdx) +{ + if (eqIdx == WQTotal) { + static const Scalar s = Parameters::Get>(); + return s; + } + if (eqIdx == SPres) { + static const Scalar s = Parameters::Get>(); + return s; + } + return Scalar{1}; // fractions and temperature stay unscaled +} + template void MultisegmentWellPrimaryVariables:: setEvaluationsFromValues() @@ -60,7 +78,10 @@ setEvaluationsFromValues() for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { evaluation_[seg][eq_idx] = 0.0; evaluation_[seg][eq_idx].setValue(value_[seg][eq_idx]); - evaluation_[seg][eq_idx].setDerivative(eq_idx + Indices::numEq, 1.0); + // 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)); } } } @@ -196,8 +217,10 @@ updateNewton(const BVectorWell& dwells, // update the segment pressure { - const int sign = dwells[seg][SPres] > 0.? 1 : -1; - const Scalar dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); + // the solution is in scaled units, the limit in physical ones + const Scalar dspres = varScale(SPres) * dwells[seg][SPres]; + const int sign = dspres > 0.? 1 : -1; + const Scalar dx_limited = sign * std::min(std::abs(dspres) * relaxation_factor, max_pressure_change); // some cases might have defaulted bhp constraint of 1 bar, we use a slightly smaller value as the bhp lower limit for Newton update // so that bhp constaint can be an active control when needed. const Scalar lower_limit = (seg == 0) ? bhp_lower_limit : seg_pres_lower_limit; @@ -206,7 +229,8 @@ updateNewton(const BVectorWell& dwells, // update the total rate // TODO: should we have a limitation of the total rate change? { - value_[seg][WQTotal] = old_primary_variables[seg][WQTotal] - relaxation_factor * dwells[seg][WQTotal]; + value_[seg][WQTotal] = old_primary_variables[seg][WQTotal] + - relaxation_factor * varScale(WQTotal) * dwells[seg][WQTotal]; // make sure that no injector produce and no producer inject if (seg == 0) { diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp index 89bff6c1e8e..0215c32cf02 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp @@ -160,6 +160,10 @@ class MultisegmentWellPrimaryVariables DeferredLogger& deferred_logger) const; private: + //! \brief Scaling of well primary variable \p eqIdx, 1.0 unless configured. + //! \details Only the Jacobian column is scaled; value_ stays physical. + static Scalar varScale(const int eqIdx); + //! \brief Initialize evaluations from values. void setEvaluationsFromValues(); From 027f004e260685e0a00060c79bb6df3441e60658 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 12 Aug 2026 09:26:03 +0200 Subject: [PATCH 3/3] Name the solver-to-physical conversion and pin the units contract 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 --- .../MultisegmentWellPrimaryVariables.cpp | 5 ++-- .../MultisegmentWellPrimaryVariables.hpp | 21 +++++++++++++- .../wells/StandardWellPrimaryVariables.cpp | 5 ++-- .../wells/StandardWellPrimaryVariables.hpp | 18 +++++++++++- tests/test_wellmodel.cpp | 29 +++++++++++++++++++ 5 files changed, 70 insertions(+), 8 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 63046daa5e1..9892f0b6fcc 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -217,8 +217,7 @@ updateNewton(const BVectorWell& dwells, // update the segment pressure { - // the solution is in scaled units, the limit in physical ones - const Scalar dspres = varScale(SPres) * dwells[seg][SPres]; + const Scalar dspres = this->physicalIncrement(dwells, seg, SPres); const int sign = dspres > 0.? 1 : -1; const Scalar dx_limited = sign * std::min(std::abs(dspres) * relaxation_factor, max_pressure_change); // some cases might have defaulted bhp constraint of 1 bar, we use a slightly smaller value as the bhp lower limit for Newton update @@ -230,7 +229,7 @@ updateNewton(const BVectorWell& dwells, // update the total rate // TODO: should we have a limitation of the total rate change? { value_[seg][WQTotal] = old_primary_variables[seg][WQTotal] - - relaxation_factor * varScale(WQTotal) * dwells[seg][WQTotal]; + - relaxation_factor * this->physicalIncrement(dwells, seg, WQTotal); // make sure that no injector produce and no producer inject if (seg == 0) { diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp index 0215c32cf02..7cfe654d7ab 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp @@ -161,9 +161,16 @@ class MultisegmentWellPrimaryVariables private: //! \brief Scaling of well primary variable \p eqIdx, 1.0 unless configured. - //! \details Only the Jacobian column is scaled; value_ stays physical. static Scalar varScale(const int eqIdx); + //! \brief The physical increment for \p eqIdx in \p seg from a solver-space update. + //! \details The only place solver-space quantities enter this class, see the + //! note on value_ below. + Scalar physicalIncrement(const BVectorWell& dwells, + const int seg, + const int eqIdx) const + { return varScale(eqIdx) * dwells[seg][eqIdx]; } + //! \brief Initialize evaluations from values. void setEvaluationsFromValues(); @@ -176,6 +183,18 @@ class MultisegmentWellPrimaryVariables //! \brief The values for the primary variables //! \details Based on different solution strategies, the wells can have different primary variables + //! + //! Units contract when scaling is enabled (varScale() != 1): + //! - value_ is PHYSICAL, and so is everything exchanged with WellState and + //! every segment pressure or rate handed to a PVT lookup. + //! - The linear system's unknown is X = x/varScale, so dwells arrives in + //! X-space; physicalIncrement() is the only conversion, and the absolute + //! limits in updateNewton() are therefore compared in physical units. + //! - eval(seg)[i] has a physical value and derivative varScale(i), so any + //! Jacobian entry must come from a derivative rather than a literal. + //! An alternative is to store X here and multiply out in the accessors; that + //! keeps value_ and the solution in one space but moves the conversion to + //! every WellState exchange and every getter. std::vector> value_; //! \brief The Evaluation for the well primary variables. diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index bc31c5e4f2b..af8a88ac74f 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp @@ -323,8 +323,7 @@ updateNewton(const BVectorWell& dwells, this->processFractions(); // updating the total rates Q_t - // The solution is in scaled units, value_ is physical, so convert the increment. - value_[WQTotal] -= varScale(WQTotal) * dwells[0][WQTotal]; + value_[WQTotal] -= this->physicalIncrement(dwells, WQTotal); // here, we make sure it is zero for wells with zero rate target(including stopped wells) if (stop_or_zero_rate_target) { @@ -339,7 +338,7 @@ updateNewton(const BVectorWell& dwells, } // updating the bottom hole pressure - const Scalar dbhp = varScale(Bhp) * dwells[0][Bhp]; + const Scalar dbhp = this->physicalIncrement(dwells, Bhp); const int sign1 = dbhp > 0 ? 1: -1; const Scalar dx1_limited = sign1 * std::min(std::abs(dbhp), std::abs(value_[Bhp]) * dBHPLimit); diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index 6364903620c..64b4c067ef8 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -155,9 +155,14 @@ class StandardWellPrimaryVariables { private: //! \brief Scaling of well primary variable \p eqIdx, 1.0 unless configured. - //! \details Only the Jacobian column is scaled; value_ stays physical. static Scalar varScale(const int eqIdx); + //! \brief The physical increment for \p eqIdx from a solver-space update. + //! \details The only place solver-space quantities enter this class, see the + //! note on value_ below. + Scalar physicalIncrement(const BVectorWell& dwells, const int eqIdx) const + { return varScale(eqIdx) * dwells[0][eqIdx]; } + //! \brief Initialize evaluations from values. void setEvaluationsFromValues(); @@ -174,6 +179,17 @@ class StandardWellPrimaryVariables { //! \brief The values for the primary variables. //! \details Based on different solution strategies, the wells can have different primary variables. + //! + //! Units contract when scaling is enabled (varScale() != 1): + //! - value_ is PHYSICAL, and so is everything exchanged with WellState. + //! - The linear system's unknown is X = x/varScale, so dwells/xw arrive in + //! X-space; physicalIncrement() is the only conversion, and the absolute + //! limits in updateNewton() are therefore compared in physical units. + //! - 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. + //! An alternative is to store X here and multiply out in the accessors; that + //! keeps value_ and the solution in one space but moves the conversion to + //! every WellState exchange and every getter. std::vector value_; //! \brief The Evaluation for the well primary variables. diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index 9dcceb2fca8..2520f49e4c5 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -252,6 +252,7 @@ BOOST_AUTO_TEST_CASE(TestPrimaryVariableScaling) { const StandardWell well(well_ecl, pinfo, setup_test.current_timestep, param, rateConverter, 0, 3, 3, 0, pdata); + constexpr double pv_bhp_scale = 65536.0; using PV = std::decay_t; PV pv(well); pv.resize(well.numStaticWellEq); @@ -279,4 +280,32 @@ BOOST_AUTO_TEST_CASE(TestPrimaryVariableScaling) { // The dimensionless fractions stay unscaled. BOOST_CHECK_EQUAL(pv.eval(PV::WFrac).derivative(numEq + PV::WFrac), 1.0); BOOST_CHECK_EQUAL(pv.eval(PV::GFrac).derivative(numEq + PV::GFrac), 1.0); + + // Every slot: the stored value and the Evaluation agree, and only the own + // column carries a derivative. This is the assertion a forgotten conversion + // in setEvaluationsFromValues would break. + for (int i = 0; i < well.numStaticWellEq; ++i) { + BOOST_CHECK_EQUAL(pv.eval(i).value(), pv.value(i)); + for (int j = 0; j < well.numStaticWellEq; ++j) { + if (i != j) { + BOOST_CHECK_EQUAL(pv.eval(i).derivative(numEq + j), 0.0); + } + } + } + + // value()/setValue() must be exact inverses, or the getPrimaryVars round + // trip used by NLDD would rescale on every save/restore. + const double bhp_phys = 300.0 * Opm::unit::barsa; + pv.setValue(PV::Bhp, bhp_phys); + BOOST_CHECK_EQUAL(pv.value(PV::Bhp), bhp_phys); + + // Newton step with the absolute bhp floor ACTIVE. The limit is physical, the + // increment arrives scaled; a missing conversion on either side lands + // somewhere other than exactly the floor. The unscaled test above never + // reaches this branch because its increment is zero. + constexpr double bhp_floor = 1.0 * Opm::unit::barsa - 1.0 * Opm::unit::Pascal; + pv.setValue(PV::Bhp, 1.5 * Opm::unit::barsa); + dwells[0][PV::Bhp] = (1.0 * Opm::unit::barsa) / pv_bhp_scale; // scaled: drives well below the floor + pv.updateNewton(dwells, false, 0.2, 1.0, logger); + BOOST_CHECK_EQUAL(pv.value(PV::Bhp), bhp_floor); }