Merge pull request #4620 from hakonhagland/wvfpdp

Add support for WVFPDP

ebos/eclgenericvanguard.cc

@@ -63,6 +63,7 @@
 #include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
 #include <opm/input/eclipse/Schedule/Well/WListManager.hpp>
+#include <opm/input/eclipse/Schedule/Well/WVFPDP.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPEXP.hpp>
 #include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 #include <opm/input/eclipse/Python/Python.hpp>

opm/simulators/utils/ParallelSerialization.cpp

@@ -57,6 +57,7 @@
 #include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
 #include <opm/input/eclipse/Schedule/Well/WList.hpp>
 #include <opm/input/eclipse/Schedule/Well/WListManager.hpp>
+#include <opm/input/eclipse/Schedule/Well/WVFPDP.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPEXP.hpp>
 
 #include <ebos/eclmpiserializer.hh>

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -754,7 +754,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"WTEMPQ", {true, std::nullopt}},
         {"WTHPMAX", {true, std::nullopt}},
         {"WTMULT", {true, std::nullopt}},
-        {"WVFPDP", {true, std::nullopt}},
         {"WWPAVE", {true, std::nullopt}},
         {"ZIPPY2", {false, std::nullopt}},
         {"ZIPP2OFF", {false, std::nullopt}},

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -209,6 +209,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
                 const double rho,
                 const bool stop_or_zero_rate_target,
                 WellState& well_state,
+                const SummaryState& summary_state,
                 DeferredLogger& deferred_logger) const
 {
     static constexpr int Gas = BlackoilPhases::Vapour;
@@ -400,7 +401,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
                    {FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
                     FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
                     FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
-                   well_state, deferred_logger);
+                   well_state, summary_state, deferred_logger);
 }
 
 template<class FluidSystem, class Indices, class Scalar>

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -25,6 +25,7 @@
 #include <opm/material/densead/Evaluation.hpp>
 
 #include <opm/simulators/wells/MultisegmentWellEquations.hpp>
+#include <opm/input/eclipse/Schedule/SummaryState.hpp>
 
 #include <array>
 #include <vector>
@@ -102,6 +103,7 @@ public:
                          const double rho,
                          const bool stop_or_zero_rate_target,
                          WellState& well_state,
+                         const SummaryState& summary_state,
                          DeferredLogger& deferred_logger) const;
 
     //! \brief Returns scaled volume fraction for a component in a segment.

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -645,7 +645,7 @@ namespace Opm
                                               max_pressure_change);
 
         this->primary_variables_.copyToWellState(*this, getRefDensity(), stop_or_zero_rate_target,
-                                                 well_state, deferred_logger);
+                                                 well_state, summary_state, deferred_logger);
         Base::calculateReservoirRates(well_state.well(this->index_of_well_));
     }
 
@@ -1135,12 +1135,13 @@ namespace Opm
                 std::vector<double> well_rates_bhp_limit;
                 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
 
+                const double thp_limit = this->getTHPConstraint(summaryState);
                 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
                                                                                    bhp_limit,
                                                                                    this->getRefDensity(),
                                                                                    this->wellEcl().alq_value(),
+                                                                                   thp_limit,
                                                                                    deferred_logger);
-                const double thp_limit = this->getTHPConstraint(summaryState);
                 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
                     this->operability_status_.obey_thp_limit_under_bhp_limit = false;
                 }

opm/simulators/wells/StandardWell.hpp

@@ -350,6 +350,7 @@ namespace Opm
 
         void updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
                                                  WellState& well_state,
+                                                 const SummaryState& summary_state,
                                                  DeferredLogger& deferred_logger) const;
 
         virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,

opm/simulators/wells/StandardWellEval.cpp

@@ -73,6 +73,7 @@ void
 StandardWellEval<FluidSystem,Indices,Scalar>::
 updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
                                     WellState& well_state,
+                                    const SummaryState& summary_state,
                                     DeferredLogger& deferred_logger) const
 {
     this->primary_variables_.copyToWellState(well_state, deferred_logger);
@@ -84,7 +85,7 @@ updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
                       {FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
                        FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
                        FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
-                      well_state, deferred_logger);
+                      well_state, summary_state, deferred_logger);
 }
 
 template<class FluidSystem, class Indices, class Scalar>

opm/simulators/wells/StandardWellEval.hpp

@@ -95,6 +95,7 @@ protected:
 
     void updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
                                              WellState& well_state,
+                                             const SummaryState& summary_state,
                                              DeferredLogger& deferred_logger) const;
 
     PrimaryVariables primary_variables_; //!< Primary variables for well

opm/simulators/wells/StandardWell_impl.hpp

@@ -917,7 +917,7 @@ namespace Opm
         const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
         updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
 
-        updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, deferred_logger);
+        updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
         Base::calculateReservoirRates(well_state.well(this->index_of_well_));
     }
 
@@ -953,9 +953,10 @@ namespace Opm
     StandardWell<TypeTag>::
     updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
                                         WellState& well_state,
+                                        const SummaryState& summary_state,
                                         DeferredLogger& deferred_logger) const
     {
-        this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, deferred_logger);
+        this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
 
         // other primary variables related to polymer injectivity study
         if constexpr (Base::has_polymermw) {
@@ -1096,12 +1097,13 @@ namespace Opm
                 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
 
                 this->adaptRatesForVFP(well_rates_bhp_limit);
+                const double thp_limit = this->getTHPConstraint(summaryState);
                 const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
                                                                                    bhp_limit,
                                                                                    this->connections_.rho(),
                                                                                    this->getALQ(well_state),
+                                                                                   thp_limit,
                                                                                    deferred_logger);
-                const double thp_limit = this->getTHPConstraint(summaryState);
                 if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
                     this->operability_status_.obey_thp_limit_under_bhp_limit = false;
                 }

opm/simulators/wells/WellBhpThpCalculator.cpp

@@ -28,8 +28,8 @@
 
 #include <opm/input/eclipse/Schedule/VFPInjTable.hpp>
 #include <opm/input/eclipse/Schedule/Well/Well.hpp>
-
 #include <opm/input/eclipse/Units/Units.hpp>
+#include <opm/input/eclipse/Schedule/Well/WVFPDP.hpp>
 
 #include <opm/material/densead/Evaluation.hpp>
 
@@ -43,7 +43,10 @@
 #include <cmath>
 #include <optional>
 
+#include <fmt/format.h>
+
 static constexpr bool extraBhpAtThpLimitOutput = false;
+static constexpr bool extraThpFromBhpOutput = false;
 
 namespace Opm
 {
@@ -103,6 +106,7 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rate
                                                  const double bhp,
                                                  const double rho,
                                                  const std::optional<double>& alq,
+                                                 const double thp_limit,
                                                  DeferredLogger& deferred_logger) const
 {
     assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
@@ -117,26 +121,70 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rate
 
     // pick the density in the top layer
     double thp = 0.0;
+    const int table_id = well_.wellEcl().vfp_table_number();
     if (well_.isInjector()) {
         assert(!alq.has_value());
-        const int table_id = well_.wellEcl().vfp_table_number();
         const double vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
         const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
-
+        auto thp_func =
-        thp = well_.vfpProperties()->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
+            [this, table_id, vfp_ref_depth, aqua, liquid, vapour, dp](
+                const double bhp_value, const double pressure_loss) {
+                    return this->well_.vfpProperties()->getInj()->thp(
+                           table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss);
+                };
+        thp = findThpFromBhpIteratively(thp_func, bhp, thp_limit, dp, deferred_logger);
     }
     else if (well_.isProducer()) {
-        const int table_id = well_.wellEcl().vfp_table_number();
         const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
         const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
-
+        auto thp_func =
-        assert(alq.has_value());
+            [this, table_id, vfp_ref_depth, aqua, liquid, vapour, dp, &alq]
-        thp = well_.vfpProperties()->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq.value());
+                (const double bhp_value, const double pressure_loss) {
+                    return this->well_.vfpProperties()->getProd()->thp(
+                       table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss, alq.value());
+                };
+        thp = findThpFromBhpIteratively(thp_func, bhp, thp_limit, dp, deferred_logger);
     }
     else {
         OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
     }
+    return thp;
+}
 
+double
+WellBhpThpCalculator::
+findThpFromBhpIteratively(
+    const std::function<double(const double, const double)>& thp_func,
+    const double bhp,
+    const double thp_limit,
+    const double dp,
+    DeferredLogger& deferred_logger) const
+{
+    auto pressure_loss = getVfpBhpAdjustment(bhp + dp, thp_limit);
+    auto thp = thp_func(bhp, pressure_loss);
+    const double tolerance = 1e-5 * unit::barsa;
+    bool do_iterate = true;
+    int it = 1;
+    int max_iterations = 50;
+    while(do_iterate) {
+        if (it > max_iterations) {
+            break;
+        }
+        double thp_prev = thp;
+        pressure_loss = getVfpBhpAdjustment(bhp + dp - pressure_loss, thp_prev);
+        thp = thp_func(bhp, pressure_loss);
+        auto error = std::fabs(thp-thp_prev);
+        if (extraThpFromBhpOutput) {
+            const std::string msg = fmt::format(
+                "findThpFromBhpIteratively(): iteration {}, thp = {}, bhp = {}, "
+                "pressure_loss = {}, error = {}", it, thp, bhp+dp-pressure_loss, pressure_loss, error);
+            deferred_logger.debug(msg);
+        }
+        if (std::fabs(thp-thp_prev) < tolerance) {
+            break;
+        }
+        it++;
+    }
     return thp;
 }
 
@@ -181,8 +229,16 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
         const auto& wfr =  well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
         const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
         const bool use_vfpexp = well_.useVfpExplicit();
-        return well_.vfpProperties()->getProd()
+        const double bhp = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
-        ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value, wfr, gfr, use_vfpexp) - dp;
+                                                     rates[Water],
+                                                     rates[Oil],
+                                                     rates[Gas],
+                                                     thp_limit,
+                                                     alq_value,
+                                                     wfr,
+                                                     gfr,
+                                                     use_vfpexp);
+        return bhp - dp + getVfpBhpAdjustment(bhp, thp_limit);
     };
 
     // Make the flo() function.
@@ -231,6 +287,7 @@ void WellBhpThpCalculator::updateThp(const double rho,
                                      const std::function<double()>& alq_value,
                                      const std::array<unsigned,3>& active,
                                      WellState& well_state,
+                                     const SummaryState& summary_state,
                                      DeferredLogger& deferred_logger) const
 {
     static constexpr int Gas = BlackoilPhases::Vapour;
@@ -263,9 +320,9 @@ void WellBhpThpCalculator::updateThp(const double rho,
     if (active[Gas]) {
         rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
     }
-
     const std::optional<double> alq = this->well_.isProducer() ? std::optional<double>(alq_value()) : std::nullopt;
-    ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, alq, deferred_logger);
+    const double thp_limit = well_.getTHPConstraint(summary_state);
+    ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, alq, thp_limit, deferred_logger);
 }
 
 template<class EvalWell>
@@ -292,32 +349,47 @@ calculateBhpFromThp(const WellState& well_state,
     const EvalWell aqua = rates[Water];
     const EvalWell liquid = rates[Oil];
     const EvalWell vapour = rates[Gas];
-
+    const double thp_limit = well_.getTHPConstraint(summaryState);
-    // pick the reference density
+    double vfp_ref_depth;
-    // typically the reference in the top layer
+    EvalWell bhp_tab;
     if (well_.isInjector() )
     {
         const auto& controls = well.injectionControls(summaryState);
-        const double vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number).getDatumDepth();
+        vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number).getDatumDepth();
-        const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
+        bhp_tab = well_.vfpProperties()->getInj()->bhp(
-        return well_.vfpProperties()->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, well_.getTHPConstraint(summaryState)) - dp;
+               controls.vfp_table_number, aqua, liquid, vapour, thp_limit);
-     }
+    }
-     else if (well_.isProducer()) {
+    else if (well_.isProducer()) {
-         const auto& controls = well.productionControls(summaryState);
+        const auto& controls = well.productionControls(summaryState);
-         const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
+        vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
-         const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
+        const auto& wfr =  well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
-         const auto& wfr =  well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
+        const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
-         const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
+        const bool use_vfpexplicit = well_.useVfpExplicit();
-         const bool use_vfpexplicit = well_.useVfpExplicit();
+        bhp_tab = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
-         return well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
                                                       aqua, liquid, vapour,
-                                                      well_.getTHPConstraint(summaryState),
+                                                      thp_limit,
                                                       well_.getALQ(well_state),
-                                                      wfr, gfr, use_vfpexplicit) - dp;
+                                                      wfr, gfr, use_vfpexplicit);
-     }
+    }
-     else {
+    else {
-         OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger);
+        OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger);
-     }
+    }
+    double bhp_tab_double_value;
+    if constexpr (std::is_same_v<EvalWell, double>) {
+        bhp_tab_double_value = bhp_tab;
+    }
+    else {  // EvalWell and bhp_tab is of type DenseAd::Evaluation<double,...,...>
+        bhp_tab_double_value = bhp_tab.value();
+    }
+    const auto bhp_adjustment = getVfpBhpAdjustment(bhp_tab_double_value, thp_limit);
+    const double dp_hydro = wellhelpers::computeHydrostaticCorrection(
+            well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
+    return bhp_tab - dp_hydro + bhp_adjustment;
+}
+
+double WellBhpThpCalculator::getVfpBhpAdjustment(const double bhp_tab, const double thp_limit) const
+{
+    return well_.wellEcl().getWVFPDP().getPressureLoss(bhp_tab, thp_limit);
 }
 
 template<class ErrorPolicy>
@@ -383,8 +455,9 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
                                                                 rho, well_.gravity());
     auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
         assert(rates.size() == 3);
-        return well_.vfpProperties()->getInj()
+        const auto bhp = well_.vfpProperties()->getInj()
-                ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp;
+                ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit);
+        return bhp - dp + getVfpBhpAdjustment(bhp, thp_limit);
     };
 
     // Make the flo() function.

opm/simulators/wells/WellBhpThpCalculator.hpp

@@ -58,6 +58,7 @@ public:
                                const double bhp,
                                const double rho,
                                const std::optional<double>& alq,
+                               const double thp_limit,
                                DeferredLogger& deferred_logger) const;
 
     //! \brief Compute BHP from THP limit for a producer.
@@ -86,6 +87,7 @@ public:
                    const std::function<double()>& alq_value,
                    const std::array<unsigned,3>& active,
                    WellState& well_state,
+                   const SummaryState& summary_state,
                    DeferredLogger& deferred_logger) const;
 
   template<class EvalWell>
@@ -122,6 +124,9 @@ private:
                          const std::array<double, 2>& range,
                          DeferredLogger& deferred_logger) const;
 
+    //! \brief Get pressure adjustment to the bhp calculated from VFP table
+    double getVfpBhpAdjustment(const double bph_tab, const double thp_limit) const;
+
     //! \brief Find limits using bisection.
     bool bisectBracket(const std::function<double(const double)>& eq,
                        const std::array<double, 2>& range,
@@ -135,6 +140,12 @@ private:
                                   double& low, double& high,
                                   DeferredLogger& deferred_logger);
 
+    double findThpFromBhpIteratively(const std::function<double(const double, const double)>& thp_func,
+                                     const double bhp,
+                                     const double thp_limit,
+                                     const double dp,
+                                     DeferredLogger& deferred_logger) const;
+
     const WellInterfaceGeneric& well_; //!< Reference to well interface
 };
 

tests/test_ParallelSerialization.cpp

@@ -108,6 +108,7 @@
 #include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
 #include <opm/input/eclipse/Schedule/Well/WList.hpp>
 #include <opm/input/eclipse/Schedule/Well/WListManager.hpp>
+#include <opm/input/eclipse/Schedule/Well/WVFPDP.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPEXP.hpp>
 #include <opm/input/eclipse/Schedule/WriteRestartFileEvents.hpp>
 #include <opm/input/eclipse/EclipseState/SimulationConfig/BCConfig.hpp>