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/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..9892f0b6fcc 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,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) { diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp index 89bff6c1e8e..7cfe654d7ab 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/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> value_; //! \brief The Evaluation for the well primary variables. diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index ed1f89853ba..af8a88ac74f 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,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. diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index 8fb7a74f5ac..64b4c067ef8 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/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 value_; //! \brief The Evaluation for the well primary variables. diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index 70ecd438920..2520f49e4c5 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -52,6 +52,7 @@ #include #include +#include #include #include @@ -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>(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); + + constexpr double pv_bhp_scale = 65536.0; + 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); + + // 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); +}