Add optional scaling of the well bhp and rate primary variables

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -133,6 +133,10 @@ void BlackoilModelParameters<Scalar>::registerParameters()
 {
     Parameters::Register<Parameters::DbhpMaxRel<Scalar>>
         ("Maximum relative change of the bottom-hole pressure in a single iteration");
+    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");
     Parameters::Register<Parameters::DwellFractionMax<Scalar>>
         ("Maximum absolute change of a well's volume fraction in a single iteration");
     Parameters::Register<Parameters::InjMultOscThreshold<Scalar>>

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -93,6 +93,19 @@ struct ToleranceWells { static constexpr Scalar value = 1e-4; };
 template<class Scalar>
 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<class Scalar>
+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<class Scalar>
+struct WellRateScaling { static constexpr Scalar value = 1.0; };
+
 struct MaxWelleqIter { static constexpr int value = 30; };
 
 template<class Scalar>

opm/simulators/wells/MultisegmentWellEquations.cpp

@@ -35,6 +35,9 @@
 #include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
 #endif
 
+#include <opm/models/utils/parametersystem.hpp>
+
+#include <opm/simulators/flow/BlackoilModelParameters.hpp>
 #include <opm/simulators/linalg/istlsparsematrixadapter.hh>
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
@@ -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<Parameters::WellBhpScaling<Scalar>>();
+        diag_ell *= bhp_scale;
+
 #define EXTRA_DEBUG_MSW 0
 #if EXTRA_DEBUG_MSW
         if (diag_ell <= 0.0) {

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -31,6 +31,9 @@
 #include <opm/models/blackoil/blackoilonephaseindices.hh>
 #include <opm/models/blackoil/blackoiltwophaseindices.hh>
 
+#include <opm/models/utils/parametersystem.hpp>
+
+#include <opm/simulators/flow/BlackoilModelParameters.hpp>
 #include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
 #include <opm/simulators/wells/RateConverter.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
@@ -52,6 +55,21 @@ resize(const int numSegments)
     evaluation_.resize(numSegments);
 }
 
+template<class FluidSystem, class Indices>
+typename FluidSystem::Scalar
+MultisegmentWellPrimaryVariables<FluidSystem,Indices>::varScale(const int eqIdx)
+{
+    if (eqIdx == WQTotal) {
+        static const Scalar s = Parameters::Get<Parameters::WellRateScaling<Scalar>>();
+        return s;
+    }
+    if (eqIdx == SPres) {
+        static const Scalar s = Parameters::Get<Parameters::WellBhpScaling<Scalar>>();
+        return s;
+    }
+    return Scalar{1}; // fractions and temperature stay unscaled
+}
+
 template<class FluidSystem, class Indices>
 void MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
 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,9 @@ 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);
+            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
             // 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 +228,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 * this->physicalIncrement(dwells, seg, WQTotal);
 
             // make sure that no injector produce and no producer inject
             if (seg == 0) {

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -160,6 +160,17 @@ class MultisegmentWellPrimaryVariables
                              DeferredLogger& deferred_logger) const;
 
 private:
+    //! \brief Scaling of well primary variable \p eqIdx, 1.0 unless configured.
+    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();
 
@@ -172,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<std::array<Scalar, numWellEq>> value_;
 
     //! \brief The Evaluation for the well primary variables.

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -25,6 +25,11 @@
 #include <opm/common/Exceptions.hpp>
 #include <opm/input/eclipse/Units/Units.hpp>
 
+#include <opm/models/utils/parametersystem.hpp>
+#include <opm/simulators/flow/BlackoilModelParameters.hpp>
+
+#include <cmath>
+
 #include <dune/common/dynvector.hh>
 #include <dune/istl/bvector.hh>
 
@@ -95,18 +100,41 @@ Scalar relaxationFactorFraction(const Scalar old_value,
 
 namespace Opm {
 
+template<class FluidSystem, class Indices>
+typename FluidSystem::Scalar
+StandardWellPrimaryVariables<FluidSystem,Indices>::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<Parameters::WellRateScaling<Scalar>>();
+        return s;
+    }
+    if (eqIdx == Bhp) {
+        static const Scalar s = Parameters::Get<Parameters::WellBhpScaling<Scalar>>();
+        return s;
+    }
+    return Scalar{1}; // the fractions are dimensionless
+}
+
 template<class FluidSystem, class Indices>
 void StandardWellPrimaryVariables<FluidSystem,Indices>::
 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,7 @@ updateNewton(const BVectorWell& dwells,
     this->processFractions();
 
     // updating the total rates Q_t
-    value_[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) {
@@ -310,8 +338,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 = 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);
     // 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.

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -154,6 +154,15 @@ class StandardWellPrimaryVariables {
                              DeferredLogger& deferred_logger) const;
 
 private:
+    //! \brief Scaling of well primary variable \p eqIdx, 1.0 unless configured.
+    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();
 
@@ -170,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<Scalar> value_;
 
     //! \brief The Evaluation for the well primary variables.

tests/test_wellmodel.cpp

@@ -52,6 +52,7 @@
 #include <opm/models/utils/start.hh>
 
 #include <opm/simulators/wells/StandardWell.hpp>
+#include <opm/simulators/utils/DeferredLogger.hpp>
 #include <opm/simulators/wells/BlackoilWellModel.hpp>
 
 #include <opm/input/eclipse/Deck/Deck.hpp>
@@ -218,3 +219,93 @@ 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<Opm::Parameters::WellBhpScaling<double>>(65536.0); // 2^16
+    Opm::Parameters::SetDefault<Opm::Parameters::WellRateScaling<double>>(0.25);   // 2^-2
+
+    const SetupTest setup_test;
+    const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep);
+    const Opm::BlackoilModelParameters<double> param;
+
+    using FluidSystem = Opm::BlackOilFluidSystem<double>;
+    // Just enough initialisation for phase-usage queries to work
+    // (same device as test_rftcontainer.cpp).
+    FluidSystem::initBegin(/*numPvtRegions=*/1);
+    using RateConverterType = Opm::RateConverter::
+        SurfaceToReservoirVoidage<FluidSystem, std::vector<int>>;
+    RateConverterType rateConverter(std::vector<int>(10, 0));
+
+    const auto& well_ecl = wells_ecl[0]; // PROD1
+    Opm::PerforationData<double> dummy;
+    std::vector<Opm::PerforationData<double>> pdata(well_ecl.getConnections().size(), dummy);
+    for (auto c = 0*pdata.size(); c < pdata.size(); ++c) {
+        pdata[c].ecl_index = c;
+    }
+    Opm::ParallelWellInfo<double> pinfo{well_ecl.name()};
+    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<decltype(well.primaryVariables())>;
+    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);
+
+    // 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);
+}