From 099322b0f0b5e7b9268a47abf7a47f90f24b1847 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 19 Feb 2024 14:00:51 +0100 Subject: [PATCH 1/4] VFPInjProperties: template Scalar type --- opm/simulators/wells/VFPInjProperties.cpp | 82 +++++++++++++---------- opm/simulators/wells/VFPInjProperties.hpp | 31 +++++---- opm/simulators/wells/VFPProperties.hpp | 10 ++- 3 files changed, 66 insertions(+), 57 deletions(-) diff --git a/opm/simulators/wells/VFPInjProperties.cpp b/opm/simulators/wells/VFPInjProperties.cpp index 8f4da40f6..2da1e3749 100644 --- a/opm/simulators/wells/VFPInjProperties.cpp +++ b/opm/simulators/wells/VFPInjProperties.cpp @@ -17,45 +17,47 @@ along with OPM. If not, see . */ - - - - #include "config.h" #include +#include + #include #include -#include - #include #include namespace Opm { -double VFPInjProperties::bhp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& thp_arg) const { +template +Scalar VFPInjProperties:: +bhp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp_arg) const +{ const VFPInjTable& table = detail::getTable(m_tables, table_id); detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg); return retval.value; } -double VFPInjProperties::thp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& bhp_arg) const { +template +Scalar VFPInjProperties:: +thp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar bhp_arg) const +{ const VFPInjTable& table = detail::getTable(m_tables, table_id); //Find interpolation variables - const double flo = detail::getFlo(table, aqua, liquid, vapour); - if (std::abs(flo) < std::numeric_limits::epsilon()) { + const Scalar flo = detail::getFlo(table, aqua, liquid, vapour); + if (std::abs(flo) < std::numeric_limits::epsilon()) { return 0.; } @@ -68,34 +70,41 @@ double VFPInjProperties::thp(int table_id, * expensive, but let us assome that nthp is small */ auto flo_i = detail::findInterpData(flo, table.getFloAxis()); - std::vector bhp_array(nthp); - for (int i=0; i bhp_array(nthp); + for (int i = 0; i < nthp; ++i) { auto thp_i = detail::findInterpData(thp_array[i], thp_array); bhp_array[i] = detail::interpolate(table, flo_i, thp_i).value; } - double retval = detail::findTHP(bhp_array, thp_array, bhp_arg); - return retval; + return detail::findTHP(bhp_array, thp_array, bhp_arg); } -const VFPInjTable& VFPInjProperties::getTable(const int table_id) const { +template +const VFPInjTable& +VFPInjProperties::getTable(const int table_id) const +{ return detail::getTable(m_tables, table_id); } -bool VFPInjProperties::hasTable(const int table_id) const { +template +bool VFPInjProperties::hasTable(const int table_id) const +{ return detail::hasTable(m_tables, table_id); } -void VFPInjProperties::addTable(const VFPInjTable& new_table) { +template +void VFPInjProperties::addTable(const VFPInjTable& new_table) +{ this->m_tables.emplace( new_table.getTableNum(), new_table ); } +template template -EvalWell VFPInjProperties::bhp(const int table_id, - const EvalWell& aqua, - const EvalWell& liquid, - const EvalWell& vapour, - const double& thp) const +EvalWell VFPInjProperties::bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const Scalar thp) const { //Get the table const VFPInjTable& table = detail::getTable(m_tables, table_id); @@ -116,12 +125,14 @@ EvalWell VFPInjProperties::bhp(const int table_id, return bhp; } +template class VFPInjProperties; + #define INSTANCE(...) \ - template __VA_ARGS__ VFPInjProperties::bhp<__VA_ARGS__>(const int, \ - const __VA_ARGS__&, \ - const __VA_ARGS__&, \ - const __VA_ARGS__&, \ - const double&) const; + template __VA_ARGS__ VFPInjProperties::bhp<__VA_ARGS__>(const int, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const double) const; INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) @@ -139,4 +150,5 @@ INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) + } //Namespace Opm diff --git a/opm/simulators/wells/VFPInjProperties.hpp b/opm/simulators/wells/VFPInjProperties.hpp index d4da6d6cb..e1c7a2a72 100644 --- a/opm/simulators/wells/VFPInjProperties.hpp +++ b/opm/simulators/wells/VFPInjProperties.hpp @@ -30,9 +30,9 @@ namespace Opm { class VFPInjTable; +template class VFPInjProperties { public: - VFPInjProperties() = default; /** * Takes *no* ownership of data. */ @@ -55,11 +55,11 @@ public: * input ADB objects. */ template - EvalWell bhp(const int table_id, + EvalWell bhp(const int table_id, const EvalWell& aqua, const EvalWell& liquid, const EvalWell& vapour, - const double& thp) const; + const Scalar thp) const; /** * Returns the table associated with the ID, or throws an exception if @@ -75,7 +75,8 @@ public: /** * Returns true if no vfp tables are in the current map */ - bool empty() const { + bool empty() const + { return m_tables.empty(); } @@ -90,11 +91,11 @@ public: * @return The bottom hole pressure, interpolated/extrapolated linearly using * the above parameters from the values in the input table. */ - double bhp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& thp) const; + Scalar bhp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp) const; /** * Linear interpolation of thp as a function of the input parameters @@ -107,19 +108,17 @@ public: * @return The tubing hole pressure, interpolated/extrapolated linearly using * the above parameters from the values in the input table. */ - double thp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& bhp) const; + Scalar thp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar bhp) const; protected: // Map which connects the table number with the table itself std::map> m_tables; }; - - } //namespace diff --git a/opm/simulators/wells/VFPProperties.hpp b/opm/simulators/wells/VFPProperties.hpp index 0e64a1bd0..4d5755fb3 100644 --- a/opm/simulators/wells/VFPProperties.hpp +++ b/opm/simulators/wells/VFPProperties.hpp @@ -26,7 +26,6 @@ #include #include -#include namespace Opm { @@ -61,7 +60,8 @@ public: /** * Returns the VFP properties for injection wells */ - const VFPInjProperties* getInj() const { + const VFPInjProperties* getInj() const + { return &m_inj; } @@ -93,13 +93,11 @@ public: } private: - VFPInjProperties m_inj; + VFPInjProperties m_inj; VFPProdProperties m_prod; const WellState& well_state_; - }; - -} //Namespace +} // namespace Opm #endif /* OPM_AUTODIFF_VFPPROPERTIES_HPP_ */ From 29d142b5e4d1c78ba90f603a0f27c4cfdbb0ed96 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 19 Feb 2024 14:00:51 +0100 Subject: [PATCH 2/4] VFPProdProperties: template Scalar type --- opm/simulators/wells/VFPProdProperties.cpp | 135 ++++++++++++--------- opm/simulators/wells/VFPProdProperties.hpp | 70 ++++++----- opm/simulators/wells/VFPProperties.hpp | 5 +- opm/simulators/wells/WellGroupHelpers.cpp | 2 +- opm/simulators/wells/WellGroupHelpers.hpp | 4 +- tests/test_vfpproperties.cpp | 8 +- 6 files changed, 119 insertions(+), 105 deletions(-) diff --git a/opm/simulators/wells/VFPProdProperties.cpp b/opm/simulators/wells/VFPProdProperties.cpp index ce19ba245..8114d5b03 100644 --- a/opm/simulators/wells/VFPProdProperties.cpp +++ b/opm/simulators/wells/VFPProdProperties.cpp @@ -31,21 +31,21 @@ namespace Opm { - - - -double VFPProdProperties::thp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& bhp_arg, - const double& alq) const { +template +Scalar VFPProdProperties:: +thp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar bhp_arg, + const Scalar alq) const +{ const VFPProdTable& table = detail::getTable(m_tables, table_id); // Find interpolation variables. - double flo = 0.0; - double wfr = 0.0; - double gfr = 0.0; + Scalar flo = 0.0; + Scalar wfr = 0.0; + Scalar gfr = 0.0; if (aqua == 0.0 && liquid == 0.0 && vapour == 0.0) { // All zero, likely at initial state. // Set FLO variable to minimum to avoid extrapolation. @@ -71,51 +71,56 @@ double VFPProdProperties::thp(int table_id, auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); auto alq_i = detail::findInterpData( alq, table.getALQAxis()); - std::vector bhp_array(nthp); - for (int i=0; i bhp_array(nthp); + for (int i = 0; i < nthp; ++i) { auto thp_i = detail::findInterpData(thp_array[i], thp_array); bhp_array[i] = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; } - double retval = detail::findTHP(bhp_array, thp_array, bhp_arg); - return retval; + return detail::findTHP(bhp_array, thp_array, bhp_arg); } - -double VFPProdProperties::bhp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& thp_arg, - const double& alq, - const double& explicit_wfr, - const double& explicit_gfr, - const bool use_expvfp) const { +template +Scalar VFPProdProperties:: +bhp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp_arg, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_expvfp) const +{ const VFPProdTable& table = detail::getTable(m_tables, table_id); detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq, explicit_wfr,explicit_gfr, use_expvfp); return retval.value; } - -const VFPProdTable& VFPProdProperties::getTable(const int table_id) const { +template +const VFPProdTable& +VFPProdProperties::getTable(const int table_id) const +{ return detail::getTable(m_tables, table_id); } -bool VFPProdProperties::hasTable(const int table_id) const { +template +bool VFPProdProperties::hasTable(const int table_id) const +{ return detail::hasTable(m_tables, table_id); } - -std::vector -VFPProdProperties:: -bhpwithflo(const std::vector& flos, +template +std::vector +VFPProdProperties:: +bhpwithflo(const std::vector& flos, const int table_id, - const double wfr, - const double gfr, - const double thp, - const double alq, - const double dp) const + const Scalar wfr, + const Scalar gfr, + const Scalar thp, + const Scalar alq, + const Scalar dp) const { // Get the table const VFPProdTable& table = detail::getTable(m_tables, table_id); @@ -124,7 +129,7 @@ bhpwithflo(const std::vector& flos, const auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); const auto alq_i = detail::findInterpData( alq, table.getALQAxis()); //assume constant - std::vector bhps(flos.size(), 0.); + std::vector bhps(flos.size(), 0.); for (std::size_t i = 0; i < flos.size(); ++i) { // Value of FLO is negative in OPM for producers, but positive in VFP table const auto flo_i = detail::findInterpData(-flos[i], table.getFloAxis()); @@ -137,13 +142,13 @@ bhpwithflo(const std::vector& flos, return bhps; } -double -VFPProdProperties:: +template +Scalar VFPProdProperties:: minimumBHP(const int table_id, - const double thp, - const double wfr, - const double gfr, - const double alq) const + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq) const { // Get the table const VFPProdTable& table = detail::getTable(m_tables, table_id); @@ -152,22 +157,24 @@ minimumBHP(const int table_id, return retval.second; } - - -void VFPProdProperties::addTable(const VFPProdTable& new_table) { +template +void VFPProdProperties::addTable(const VFPProdTable& new_table) +{ this->m_tables.emplace( new_table.getTableNum(), new_table ); } +template template -EvalWell VFPProdProperties::bhp(const int table_id, - const EvalWell& aqua, - const EvalWell& liquid, - const EvalWell& vapour, - const double& thp, - const double& alq, - const double& explicit_wfr, - const double& explicit_gfr, - const bool use_expvfp) const +EvalWell VFPProdProperties:: +bhp(const int table_id, + const EvalWell& aqua, + const EvalWell& liquid, + const EvalWell& vapour, + const Scalar thp, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_expvfp) const { //Get the table const VFPProdTable& table = detail::getTable(m_tables, table_id); @@ -197,10 +204,18 @@ EvalWell VFPProdProperties::bhp(const int table_id, return bhp; } +template class VFPProdProperties; + #define INSTANCE(...) \ - template __VA_ARGS__ VFPProdProperties::bhp<__VA_ARGS__>(const int, \ - const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&, \ - const double&, const double&, const double&, const double&, const bool) const; + template __VA_ARGS__ VFPProdProperties::bhp<__VA_ARGS__>(const int, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const double, \ + const double, \ + const double, \ + const double, \ + const bool) const; INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) diff --git a/opm/simulators/wells/VFPProdProperties.hpp b/opm/simulators/wells/VFPProdProperties.hpp index d5a72975a..65f6d0583 100644 --- a/opm/simulators/wells/VFPProdProperties.hpp +++ b/opm/simulators/wells/VFPProdProperties.hpp @@ -34,9 +34,9 @@ class VFPProdTable; * water fraction, gas fraction, and artificial lift for production VFP tables, and similarly * the BHP as a function of the rate and tubing head pressure. */ +template class VFPProdProperties { public: - VFPProdProperties() = default; /** * Takes *no* ownership of data. */ @@ -60,15 +60,15 @@ public: * input ADB objects. */ template - EvalWell bhp(const int table_id, + EvalWell bhp(const int table_id, const EvalWell& aqua, const EvalWell& liquid, const EvalWell& vapour, - const double& thp, - const double& alq, - const double& explicit_wfr, - const double& explicit_gfr, - const bool use_expvfp) const; + const Scalar thp, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_expvfp) const; /** * Linear interpolation of bhp as a function of the input parameters @@ -82,15 +82,15 @@ public: * @return The bottom hole pressure, interpolated/extrapolated linearly using * the above parameters from the values in the input table. */ - double bhp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& thp, - const double& alq, - const double& explicit_wfr, - const double& explicit_gfr, - const bool use_expvfp) const; + Scalar bhp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_expvfp) const; /** * Linear interpolation of thp as a function of the input parameters @@ -104,12 +104,12 @@ public: * @return The tubing hole pressure, interpolated/extrapolated linearly using * the above parameters from the values in the input table. */ - double thp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& bhp, - const double& alq) const; + Scalar thp(const int table_id, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar bhp, + const Scalar alq) const; /** * Returns the table associated with the ID, or throws an exception if @@ -125,33 +125,31 @@ public: /** * Returns true if no vfp tables are in the current map */ - bool empty() const { + bool empty() const + { return m_tables.empty(); } /** * Returns minimum bhp for given thp, wfr, gfr and alq */ - double minimumBHP(const int table_id, const double thp, - const double wfr, const double gfr, const double alq) const; + Scalar minimumBHP(const int table_id, const Scalar thp, + const Scalar wfr, const Scalar gfr, const Scalar alq) const; + protected: // calculate a group bhp values with a group of flo rate values - std::vector bhpwithflo(const std::vector& flos, + std::vector bhpwithflo(const std::vector& flos, const int table_id, - const double wfr, - const double gfr, - const double thp, - const double alq, - const double dp) const; + const Scalar wfr, + const Scalar gfr, + const Scalar thp, + const Scalar alq, + const Scalar dp) const; // Map which connects the table number with the table itself std::map> m_tables; }; - - - -} //namespace - +} // namespace Opm #endif /* OPM_AUTODIFF_VFPPRODPROPERTIES_HPP_ */ diff --git a/opm/simulators/wells/VFPProperties.hpp b/opm/simulators/wells/VFPProperties.hpp index 4d5755fb3..b6e51d0cc 100644 --- a/opm/simulators/wells/VFPProperties.hpp +++ b/opm/simulators/wells/VFPProperties.hpp @@ -68,7 +68,8 @@ public: /** * Returns the VFP properties for production wells */ - const VFPProdProperties* getProd() const { + const VFPProdProperties* getProd() const + { return &m_prod; } @@ -94,7 +95,7 @@ public: private: VFPInjProperties m_inj; - VFPProdProperties m_prod; + VFPProdProperties m_prod; const WellState& well_state_; }; diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index b2f970dd3..bb1e55856 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -811,7 +811,7 @@ WellGroupHelpers:: computeNetworkPressures(const Network::ExtNetwork& network, const WellState& well_state, const GroupState& group_state, - const VFPProdProperties& vfp_prod_props, + const VFPProdProperties& vfp_prod_props, const Schedule& schedule, const int report_time_step) { diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 47b291f1e..70754640b 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -37,7 +37,7 @@ template class GroupState; namespace Network { class ExtNetwork; } struct PhaseUsage; class Schedule; -class VFPProdProperties; +template class VFPProdProperties; template class WellState; class FieldPropsManager; @@ -202,7 +202,7 @@ public: computeNetworkPressures(const Network::ExtNetwork& network, const WellState& well_state, const GroupState& group_state, - const VFPProdProperties& vfp_prod_props, + const VFPProdProperties& vfp_prod_props, const Schedule& schedule, const int report_time_step); diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 4f33b7403..44e17b047 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -208,7 +208,7 @@ struct TrivialFixture { alq_axis, data)); - properties = std::make_shared(); + properties = std::make_shared>(); properties->addTable( *table ); } @@ -219,7 +219,7 @@ struct TrivialFixture { - std::shared_ptr properties; + std::shared_ptr> properties; std::shared_ptr table; std::vector table_ids; @@ -612,7 +612,7 @@ VFPPROD \n\ auto deck = parser.parseString(table_str); bool gaslift_active = false; Opm::VFPProdTable table(deck["VFPPROD"].front(), gaslift_active, units); - Opm::VFPProdProperties properties; + Opm::VFPProdProperties properties; properties.addTable( table ); const int n = 5; //Number of points to check per axis @@ -673,7 +673,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD) bool gaslift_active = false; Opm::VFPProdTable table(deck["VFPPROD"].front(), gaslift_active, units); - Opm::VFPProdProperties properties; + Opm::VFPProdProperties properties; properties.addTable(table); //Do some rudimentary testing From 374798134785a3b60c6906a801a4b9d1907771fd Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 19 Feb 2024 14:00:51 +0100 Subject: [PATCH 3/4] VFPProperties: template Scalar type --- examples/printvfp.cpp | 6 ++++-- .../wells/BlackoilWellModelGeneric.hpp | 4 ++-- .../wells/BlackoilWellModel_impl.hpp | 2 +- opm/simulators/wells/VFPProperties.hpp | 21 +++++++++++-------- opm/simulators/wells/WellInterfaceGeneric.cpp | 2 +- opm/simulators/wells/WellInterfaceGeneric.hpp | 8 +++---- 6 files changed, 24 insertions(+), 19 deletions(-) diff --git a/examples/printvfp.cpp b/examples/printvfp.cpp index 759258aed..a58fd6b3b 100644 --- a/examples/printvfp.cpp +++ b/examples/printvfp.cpp @@ -46,7 +46,7 @@ struct Setup std::shared_ptr python; std::unique_ptr schedule; std::unique_ptr summary_state; - std::unique_ptr vfp_properties; + std::unique_ptr> vfp_properties; Setup(const std::string& file) { @@ -63,7 +63,9 @@ struct Setup const int step = 0; const auto& sched_state = schedule->operator[](step); WellState well_state(phaseUsage(runspec.phases())); - vfp_properties = std::make_unique(sched_state.vfpinj(), sched_state.vfpprod(), well_state); + vfp_properties = std::make_unique>(sched_state.vfpinj(), + sched_state.vfpprod(), + well_state); }; }; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index 45875a5a7..e1e13c456 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -63,7 +63,7 @@ namespace Opm { class Schedule; struct SimulatorUpdate; class SummaryConfig; - class VFPProperties; + template class VFPProperties; template class WellInterfaceGeneric; template class WellState; } // namespace Opm @@ -557,7 +557,7 @@ protected: mutable std::unordered_set closed_this_step_; GuideRate guideRate_; - std::unique_ptr vfp_properties_{}; + std::unique_ptr> vfp_properties_{}; std::map node_pressures_; // Storing network pressures for output. // previous injection multiplier, it is used in the injection multiplier calculation for WINJMULT keyword diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 9044f75e6..fcfe6b258 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -275,7 +275,7 @@ namespace Opm { { const auto& sched_state = this->schedule()[timeStepIdx]; - this->vfp_properties_ = std::make_unique + this->vfp_properties_ = std::make_unique> (sched_state.vfpinj(), sched_state.vfpprod(), this->wellState()); } } diff --git a/opm/simulators/wells/VFPProperties.hpp b/opm/simulators/wells/VFPProperties.hpp index b6e51d0cc..7c41c5f1d 100644 --- a/opm/simulators/wells/VFPProperties.hpp +++ b/opm/simulators/wells/VFPProperties.hpp @@ -36,6 +36,7 @@ class VFPProdTable; * A thin wrapper class that holds one VFPProdProperties and one * VFPInjProperties object. */ +template class VFPProperties { public: /** @@ -47,8 +48,8 @@ public: VFPProperties(const std::vector>& inj_tables, const std::vector>& prod_tables, - const WellState& well_state) - :well_state_(well_state) + const WellState& well_state) + : well_state_(well_state) { for (const auto& vfpinj : inj_tables) this->m_inj.addTable( vfpinj ); @@ -60,7 +61,7 @@ public: /** * Returns the VFP properties for injection wells */ - const VFPInjProperties* getInj() const + const VFPInjProperties* getInj() const { return &m_inj; } @@ -68,12 +69,13 @@ public: /** * Returns the VFP properties for production wells */ - const VFPProdProperties* getProd() const + const VFPProdProperties* getProd() const { return &m_prod; } - double getExplicitWFR(const int table_id, const std::size_t well_index) const { + Scalar getExplicitWFR(const int table_id, const std::size_t well_index) const + { const auto& rates = well_state_.well(well_index).prev_surface_rates; const auto& pu = well_state_.phaseUsage(); const auto& aqua = pu.phase_used[BlackoilPhases::Aqua]? rates[pu.phase_pos[BlackoilPhases::Aqua]]:0.0; @@ -83,7 +85,8 @@ public: return detail::getWFR(table, aqua, liquid, vapour); } - double getExplicitGFR(const int table_id, const std::size_t well_index) const { + Scalar getExplicitGFR(const int table_id, const std::size_t well_index) const + { const auto& rates = well_state_.well(well_index).prev_surface_rates; const auto& pu = well_state_.phaseUsage(); const auto& aqua = pu.phase_used[BlackoilPhases::Aqua]? rates[pu.phase_pos[BlackoilPhases::Aqua]]:0.0; @@ -94,9 +97,9 @@ public: } private: - VFPInjProperties m_inj; - VFPProdProperties m_prod; - const WellState& well_state_; + VFPInjProperties m_inj; + VFPProdProperties m_prod; + const WellState& well_state_; }; } // namespace Opm diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 2cf2faec3..1b186bdca 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -374,7 +374,7 @@ closeCompletions(const WellTestState& wellTestState) template void WellInterfaceGeneric:: -setVFPProperties(const VFPProperties* vfp_properties_arg) +setVFPProperties(const VFPProperties* vfp_properties_arg) { vfp_properties_ = vfp_properties_arg; } diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index fa1d1eccc..715ef8290 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -40,7 +40,7 @@ class ParallelWellInfo; struct PerforationData; struct PhaseUsage; class SummaryState; -class VFPProperties; +template class VFPProperties; class WellTestState; template class WellState; template class SingleWellState; @@ -94,7 +94,7 @@ public: void initCompletions(); void closeCompletions(const WellTestState& wellTestState); - void setVFPProperties(const VFPProperties* vfp_properties_arg); + void setVFPProperties(const VFPProperties* vfp_properties_arg); void setPrevSurfaceRates(WellState& well_state, const WellState& prev_well_state) const; void setGuideRate(const GuideRate* guide_rate_arg); @@ -129,7 +129,7 @@ public: Scalar gravity() const { return gravity_; } - const VFPProperties* vfpProperties() const { return vfp_properties_; } + const VFPProperties* vfpProperties() const { return vfp_properties_; } const ParallelWellInfo& parallelWellInfo() const { return parallel_well_info_; } @@ -358,7 +358,7 @@ protected: std::vector inj_fc_multiplier_; Scalar well_efficiency_factor_; - const VFPProperties* vfp_properties_; + const VFPProperties* vfp_properties_; const GuideRate* guide_rate_; std::vector well_control_log_; From d5d16eaee4171cb52df0854e66eda84e10d6406a Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Feb 2024 12:51:56 +0100 Subject: [PATCH 4/4] VFPHelpers: move some functions into a class with static members and template Scalar type --- examples/printvfp.cpp | 12 +- opm/simulators/wells/VFPHelpers.cpp | 362 +++++++++--------- opm/simulators/wells/VFPHelpers.hpp | 218 +++++------ opm/simulators/wells/VFPInjProperties.cpp | 18 +- opm/simulators/wells/VFPProdProperties.cpp | 45 ++- opm/simulators/wells/WellBhpThpCalculator.cpp | 9 +- tests/test_vfpproperties.cpp | 14 +- 7 files changed, 355 insertions(+), 323 deletions(-) diff --git a/examples/printvfp.cpp b/examples/printvfp.cpp index a58fd6b3b..698b01478 100644 --- a/examples/printvfp.cpp +++ b/examples/printvfp.cpp @@ -83,13 +83,13 @@ double computeBhp(const VFPProdTable& table, // First, find the values to interpolate between. // Assuming positive flo here! assert(flo > 0.0); - auto flo_i = detail::findInterpData(flo, table.getFloAxis()); - auto thp_i = detail::findInterpData(thp, table.getTHPAxis()); // assume constant - auto wfr_i = detail::findInterpData(wfr, table.getWFRAxis()); - auto gfr_i = detail::findInterpData(gfr, table.getGFRAxis()); - auto alq_i = detail::findInterpData(alq, table.getALQAxis()); //assume constant + auto flo_i = VFPHelpers::findInterpData(flo, table.getFloAxis()); + auto thp_i = VFPHelpers::findInterpData(thp, table.getTHPAxis()); // assume constant + auto wfr_i = VFPHelpers::findInterpData(wfr, table.getWFRAxis()); + auto gfr_i = VFPHelpers::findInterpData(gfr, table.getGFRAxis()); + auto alq_i = VFPHelpers::findInterpData(alq, table.getALQAxis()); //assume constant - return detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; + return VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; } diff --git a/opm/simulators/wells/VFPHelpers.cpp b/opm/simulators/wells/VFPHelpers.cpp index c7c847263..ccb859b33 100644 --- a/opm/simulators/wells/VFPHelpers.cpp +++ b/opm/simulators/wells/VFPHelpers.cpp @@ -38,14 +38,15 @@ namespace { * Helper function that finds x for a given value of y for a line * *NOTE ORDER OF ARGUMENTS* */ -double findX(const double x0, - const double x1, - const double y0, - const double y1, - const double y) +template +Scalar findX(const Scalar x0, + const Scalar x1, + const Scalar y0, + const Scalar y1, + const Scalar y) { - const double dx = x1 - x0; - const double dy = y1 - y0; + const Scalar dx = x1 - x0; + const Scalar dy = y1 - y0; /** * y = y0 + (dy / dx) * (x - x0) @@ -54,7 +55,7 @@ double findX(const double x0, * If dy is zero, use x1 as the value. */ - double x = 0.0; + Scalar x = 0.0; if (dy != 0.0) { x = x0 + (y-y0) * (dx/dy); @@ -77,17 +78,18 @@ static T chopNegativeValues(const T& value) { } namespace Opm { -namespace detail { -InterpData findInterpData(const double value_in, const std::vector& values) +template +detail::InterpData VFPHelpers::findInterpData(const Scalar value_in, + const std::vector& values) { - InterpData retval; + detail::InterpData retval; const int nvalues = values.size(); // chopping the value to be zero, which means we do not // extrapolate the table towards nagative ranges - const double value = value_in < 0.? 0. : value_in; + const Scalar value = value_in < 0.? 0. : value_in; //If we only have one value in our vector, return that if (values.size() == 1) { @@ -119,8 +121,8 @@ InterpData findInterpData(const double value_in, const std::vector& valu } } - const double start = values[retval.ind_[0]]; - const double end = values[retval.ind_[1]]; + const Scalar start = values[retval.ind_[0]]; + const Scalar end = values[retval.ind_[1]]; //Find interpolation ratio if (end > start) { @@ -138,50 +140,17 @@ InterpData findInterpData(const double value_in, const std::vector& valu return retval; } -VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs) -{ - lhs.value += rhs.value; - lhs.dthp += rhs.dthp; - lhs.dwfr += rhs.dwfr; - lhs.dgfr += rhs.dgfr; - lhs.dalq += rhs.dalq; - lhs.dflo += rhs.dflo; - return lhs; -} - -VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs) -{ - lhs.value -= rhs.value; - lhs.dthp -= rhs.dthp; - lhs.dwfr -= rhs.dwfr; - lhs.dgfr -= rhs.dgfr; - lhs.dalq -= rhs.dalq; - lhs.dflo -= rhs.dflo; - return lhs; -} - -VFPEvaluation operator*(double lhs, const VFPEvaluation& rhs) -{ - VFPEvaluation retval; - retval.value = rhs.value * lhs; - retval.dthp = rhs.dthp * lhs; - retval.dwfr = rhs.dwfr * lhs; - retval.dgfr = rhs.dgfr * lhs; - retval.dalq = rhs.dalq * lhs; - retval.dflo = rhs.dflo * lhs; - return retval; -} - -VFPEvaluation interpolate(const VFPProdTable& table, - const InterpData& flo_i, - const InterpData& thp_i, - const InterpData& wfr_i, - const InterpData& gfr_i, - const InterpData& alq_i) +template +detail::VFPEvaluation VFPHelpers:: +interpolate(const VFPProdTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i, + const detail::InterpData& wfr_i, + const detail::InterpData& gfr_i, + const detail::InterpData& alq_i) { //Values and derivatives in a 5D hypercube - VFPEvaluation nn[2][2][2][2][2]; - + detail::VFPEvaluation nn[2][2][2][2][2]; //Pick out nearest neighbors (nn) to our evaluation point //This is not really required, but performance-wise it may pay off, since the 32-elements @@ -231,7 +200,7 @@ VFPEvaluation interpolate(const VFPProdTable& table, } } - double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. + Scalar t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. // Remove dimensions one by one // Example: going from 3D to 2D to 1D, we start by interpolating along @@ -280,13 +249,14 @@ VFPEvaluation interpolate(const VFPProdTable& table, return nn[0][0][0][0][0]; } -VFPEvaluation interpolate(const VFPInjTable& table, - const InterpData& flo_i, - const InterpData& thp_i) +template +detail::VFPEvaluation VFPHelpers:: +interpolate(const VFPInjTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i) { //Values and derivatives in a 2D plane - VFPEvaluation nn[2][2]; - + detail::VFPEvaluation nn[2][2]; //Pick out nearest neighbors (nn) to our evaluation point //The following ladder of for loops will presumably be unrolled by a reasonable compiler. @@ -318,7 +288,7 @@ VFPEvaluation interpolate(const VFPInjTable& table, nn[i][1].dflo = nn[i][0].dflo; } - double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. + Scalar t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. // Go from 2D to 1D t2 = flo_i.factor_; @@ -334,20 +304,22 @@ VFPEvaluation interpolate(const VFPInjTable& table, return nn[0][0]; } -VFPEvaluation bhp(const VFPProdTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp, - const double alq, - const double explicit_wfr, - const double explicit_gfr, - const bool use_vfpexplicit) +template +detail::VFPEvaluation VFPHelpers:: +bhp(const VFPProdTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_vfpexplicit) { //Find interpolation variables - double flo = detail::getFlo(table, aqua, liquid, vapour); - double wfr = detail::getWFR(table, aqua, liquid, vapour); - double gfr = detail::getGFR(table, aqua, liquid, vapour); + Scalar flo = detail::getFlo(table, aqua, liquid, vapour); + Scalar wfr = detail::getWFR(table, aqua, liquid, vapour); + Scalar gfr = detail::getGFR(table, aqua, liquid, vapour); if (use_vfpexplicit || -flo < table.getFloAxis().front()) { wfr = explicit_wfr; gfr = explicit_gfr; @@ -355,43 +327,47 @@ VFPEvaluation bhp(const VFPProdTable& table, //First, find the values to interpolate between //Recall that flo is negative in Opm, so switch sign. - auto flo_i = detail::findInterpData(-flo, table.getFloAxis()); - auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); - auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); - auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); - auto alq_i = detail::findInterpData( alq, table.getALQAxis()); + auto flo_i = findInterpData(-flo, table.getFloAxis()); + auto thp_i = findInterpData( thp, table.getTHPAxis()); + auto wfr_i = findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = findInterpData( gfr, table.getGFRAxis()); + auto alq_i = findInterpData( alq, table.getALQAxis()); - detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + detail::VFPEvaluation retval = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); return retval; } -VFPEvaluation bhp(const VFPInjTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp) +template +detail::VFPEvaluation VFPHelpers:: +bhp(const VFPInjTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp) { //Find interpolation variables - double flo = detail::getFlo(table, aqua, liquid, vapour); + Scalar flo = detail::getFlo(table, aqua, liquid, vapour); //First, find the values to interpolate between - auto flo_i = detail::findInterpData(flo, table.getFloAxis()); - auto thp_i = detail::findInterpData(thp, table.getTHPAxis()); + auto flo_i = findInterpData(flo, table.getFloAxis()); + auto thp_i = findInterpData(thp, table.getTHPAxis()); //Then perform the interpolation itself - detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i); + detail::VFPEvaluation retval = interpolate(table, flo_i, thp_i); return retval; } -double findTHP(const std::vector& bhp_array, - const std::vector& thp_array, - double bhp) +template +Scalar VFPHelpers:: +findTHP(const std::vector& bhp_array, + const std::vector& thp_array, + Scalar bhp) { int nthp = thp_array.size(); - double thp = -1e100; + Scalar thp = -1e100; //Check that our thp axis is sorted assert(std::is_sorted(thp_array.begin(), thp_array.end())); @@ -406,19 +382,19 @@ double findTHP(const std::vector& bhp_array, //Target bhp less than all values in array, extrapolate if (bhp <= bhp_array[0]) { //TODO: LOG extrapolation - const double& x0 = thp_array[0]; - const double& x1 = thp_array[1]; - const double& y0 = bhp_array[0]; - const double& y1 = bhp_array[1]; + const Scalar& x0 = thp_array[0]; + const Scalar& x1 = thp_array[1]; + const Scalar& y0 = bhp_array[0]; + const Scalar& y1 = bhp_array[1]; thp = findX(x0, x1, y0, y1, bhp); } //Target bhp greater than all values in array, extrapolate else if (bhp > bhp_array[nthp-1]) { //TODO: LOG extrapolation - const double& x0 = thp_array[nthp-2]; - const double& x1 = thp_array[nthp-1]; - const double& y0 = bhp_array[nthp-2]; - const double& y1 = bhp_array[nthp-1]; + const Scalar& x0 = thp_array[nthp-2]; + const Scalar& x1 = thp_array[nthp-1]; + const Scalar& y0 = bhp_array[nthp-2]; + const Scalar& y1 = bhp_array[nthp-1]; thp = findX(x0, x1, y0, y1, bhp); } //Target bhp within table ranges, interpolate @@ -432,8 +408,8 @@ double findTHP(const std::vector& bhp_array, int i=0; bool found = false; for (; i& bhp_array, assert(found == true); static_cast(found); //Silence compiler warning - const double& x0 = thp_array[i ]; - const double& x1 = thp_array[i+1]; - const double& y0 = bhp_array[i ]; - const double& y1 = bhp_array[i+1]; + const Scalar& x0 = thp_array[i ]; + const Scalar& x1 = thp_array[i+1]; + const Scalar& y0 = bhp_array[i ]; + const Scalar& y1 = bhp_array[i+1]; thp = findX(x0, x1, y0, y1, bhp); } } @@ -459,8 +435,8 @@ double findTHP(const std::vector& bhp_array, int i=0; bool found = false; for (; i& bhp_array, } } if (found) { - const double& x0 = thp_array[i ]; - const double& x1 = thp_array[i+1]; - const double& y0 = bhp_array[i ]; - const double& y1 = bhp_array[i+1]; + const Scalar& x0 = thp_array[i ]; + const Scalar& x1 = thp_array[i+1]; + const Scalar& y0 = bhp_array[i ]; + const Scalar& y1 = bhp_array[i+1]; thp = findX(x0, x1, y0, y1, bhp); } else if (bhp <= bhp_array[0]) { //TODO: LOG extrapolation - const double& x0 = thp_array[0]; - const double& x1 = thp_array[1]; - const double& y0 = bhp_array[0]; - const double& y1 = bhp_array[1]; + const Scalar& x0 = thp_array[0]; + const Scalar& x1 = thp_array[1]; + const Scalar& y0 = bhp_array[0]; + const Scalar& y1 = bhp_array[1]; thp = findX(x0, x1, y0, y1, bhp); } //Target bhp greater than all values in array, extrapolate else if (bhp > bhp_array[nthp-1]) { //TODO: LOG extrapolation - const double& x0 = thp_array[nthp-2]; - const double& x1 = thp_array[nthp-1]; - const double& y0 = bhp_array[nthp-2]; - const double& y1 = bhp_array[nthp-1]; + const Scalar& x0 = thp_array[nthp-2]; + const Scalar& x1 = thp_array[nthp-1]; + const Scalar& y0 = bhp_array[nthp-2]; + const Scalar& y1 = bhp_array[nthp-1]; thp = findX(x0, x1, y0, y1, bhp); } else { @@ -499,86 +475,88 @@ double findTHP(const std::vector& bhp_array, return thp; } -std::pair +template +std::pair VFPHelpers:: getMinimumBHPCoordinate(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq) -{ + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq) +{ // Given fixed thp, wfr, gfr and alq, this function finds the minimum bhp and returns - // the corresponding pair (-flo_at_bhp_min, bhp_min). No assumption is taken on the - // shape of the function bhp(flo), so all points in the flo-axis is checked. - double flo_at_bhp_min = 0.0; // start by checking flo=0 - auto flo_i = detail::findInterpData(flo_at_bhp_min, table.getFloAxis()); - auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); - auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); - auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); - auto alq_i = detail::findInterpData( alq, table.getALQAxis()); + // the corresponding pair (-flo_at_bhp_min, bhp_min). No assumption is taken on the + // shape of the function bhp(flo), so all points in the flo-axis is checked. + Scalar flo_at_bhp_min = 0.0; // start by checking flo=0 + auto flo_i = findInterpData(flo_at_bhp_min, table.getFloAxis()); + auto thp_i = findInterpData( thp, table.getTHPAxis()); + auto wfr_i = findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = findInterpData( gfr, table.getGFRAxis()); + auto alq_i = findInterpData( alq, table.getALQAxis()); - detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); - double bhp_min = bhp_i.value; + detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + Scalar bhp_min = bhp_i.value; const std::vector& flos = table.getFloAxis(); for (size_t i = 0; i < flos.size(); ++i) { - flo_i = detail::findInterpData(flos[i], flos); - bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + flo_i = findInterpData(flos[i], flos); + bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); if (bhp_i.value < bhp_min){ bhp_min = bhp_i.value; flo_at_bhp_min = flos[i]; } } // return negative flo - return std::make_pair(-flo_at_bhp_min, bhp_min); -} + return std::make_pair(-flo_at_bhp_min, bhp_min); +} -std::optional> +template +std::optional> VFPHelpers:: intersectWithIPR(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq, - const double ipr_a, - const double ipr_b, - const std::function& adjust_bhp) -{ - // Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection - // between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq, + const Scalar ipr_a, + const Scalar ipr_b, + const std::function& adjust_bhp) +{ + // Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection + // between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable // intersections, the one corresponding the largest flo is returned. // The adjust_bhp-function is used to adjust the vfp-table bhp-values to actual bhp-values due - // vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments. + // vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments. // NOTE: ipr-line is q=b*bhp - a! // ipr is given for negative flo, so - // flo = -b*bhp + a, i.e., bhp = -(flo-a)/b - auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); - auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); - auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); - auto alq_i = detail::findInterpData( alq, table.getALQAxis()); + // flo = -b*bhp + a, i.e., bhp = -(flo-a)/b + auto thp_i = findInterpData( thp, table.getTHPAxis()); + auto wfr_i = findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = findInterpData( gfr, table.getGFRAxis()); + auto alq_i = findInterpData( alq, table.getALQAxis()); if (ipr_b == 0.0) { // this shouldn't happen, but deal with it to be safe - auto flo_i = detail::findInterpData(ipr_a, table.getFloAxis()); - detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + auto flo_i = findInterpData(ipr_a, table.getFloAxis()); + detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); return std::make_pair(-ipr_a, adjust_bhp(bhp_i.value)); } // find largest flo (flo_x) for which y = bhp(flo) + (flo-a)/b = 0 and dy/dflo > 0 - double flo_x = -1.0; - double flo0, flo1; - double y0, y1; + Scalar flo_x = -1.0; + Scalar flo0, flo1; + Scalar y0, y1; flo0 = 0.0; // start by checking flo=0 - auto flo_i = detail::findInterpData(flo0, table.getFloAxis()); - detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + auto flo_i = findInterpData(flo0, table.getFloAxis()); + detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); y0 = adjust_bhp(bhp_i.value) - ipr_a/ipr_b; // +0.0/ipr_b - std::vector flos = table.getFloAxis(); + const std::vector& flos = table.getFloAxis(); for (size_t i = 0; i < flos.size(); ++i) { flo1 = flos[i]; - flo_i = detail::findInterpData(flo1, table.getFloAxis()); - bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + flo_i = findInterpData(flo1, flos); + bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); y1 = adjust_bhp(bhp_i.value) + (flo1 - ipr_a)/ipr_b; if (y0 < 0 && y1 >= 0){ // crossing with positive slope - double w = -y0/(y1-y0); + Scalar w = -y0/(y1-y0); w = std::clamp(w, 0.0, 1.0); // just to be safe (if y0~y1~0) flo_x = flo0 + w*(flo1 - flo0); } @@ -591,7 +569,46 @@ intersectWithIPR(const VFPProdTable& table, } else { return std::nullopt; } -} +} + +namespace detail { + +template +VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs) +{ + lhs.value += rhs.value; + lhs.dthp += rhs.dthp; + lhs.dwfr += rhs.dwfr; + lhs.dgfr += rhs.dgfr; + lhs.dalq += rhs.dalq; + lhs.dflo += rhs.dflo; + return lhs; +} + +template +VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs) +{ + lhs.value -= rhs.value; + lhs.dthp -= rhs.dthp; + lhs.dwfr -= rhs.dwfr; + lhs.dgfr -= rhs.dgfr; + lhs.dalq -= rhs.dalq; + lhs.dflo -= rhs.dflo; + return lhs; +} + +template +VFPEvaluation operator*(Scalar lhs, const VFPEvaluation& rhs) +{ + VFPEvaluation retval; + retval.value = rhs.value * lhs; + retval.dthp = rhs.dthp * lhs; + retval.dwfr = rhs.dwfr * lhs; + retval.dgfr = rhs.dgfr * lhs; + retval.dalq = rhs.dalq * lhs; + retval.dflo = rhs.dflo * lhs; + return retval; +} template T getFlo(const VFPProdTable& table, @@ -753,4 +770,7 @@ INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) } // namespace detail + +template class VFPHelpers; + } // namespace Opm diff --git a/opm/simulators/wells/VFPHelpers.hpp b/opm/simulators/wells/VFPHelpers.hpp index df678da98..761035ee5 100644 --- a/opm/simulators/wells/VFPHelpers.hpp +++ b/opm/simulators/wells/VFPHelpers.hpp @@ -21,7 +21,6 @@ #ifndef OPM_AUTODIFF_VFPHELPERS_HPP_ #define OPM_AUTODIFF_VFPHELPERS_HPP_ -#include #include #include #include @@ -37,6 +36,40 @@ class VFPProdTable; namespace detail { +/** + * An "ADB-like" structure with a single value and a set of derivatives + */ +template +struct VFPEvaluation +{ + VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {}; + Scalar value; + Scalar dthp; + Scalar dwfr; + Scalar dgfr; + Scalar dalq; + Scalar dflo; +}; + +template +VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs); +template +VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs); +template +VFPEvaluation operator*(Scalar lhs, const VFPEvaluation& rhs); + +/** + * Helper struct for linear interpolation + */ +template +struct InterpData +{ + InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {} + int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value] + Scalar inv_dist_; // 1 / distance between the two end points of the segment. Used to calculate derivatives and uses 1.0 / 0.0 = 0.0 as a convention + Scalar factor_; // Interpolation factor +}; + /** * Computes the flo parameter according to the flo_type_ * for production tables @@ -79,76 +112,6 @@ T getGFR(const VFPProdTable& table, const T& liquid, const T& vapour); -/** - * Helper struct for linear interpolation - */ -struct InterpData { - InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {} - int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value] - double inv_dist_; // 1 / distance between the two end points of the segment. Used to calculate derivatives and uses 1.0 / 0.0 = 0.0 as a convention - double factor_; // Interpolation factor -}; - -/** - * Helper function to find indices etc. for linear interpolation and extrapolation - * @param value_in Value to find in values - * @param values Sorted list of values to search for value in. - * @return Data required to find the interpolated value - */ -InterpData findInterpData(const double value_in, const std::vector& values); - -/** - * An "ADB-like" structure with a single value and a set of derivatives - */ -struct VFPEvaluation { - VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {}; - double value; - double dthp; - double dwfr; - double dgfr; - double dalq; - double dflo; -}; - -VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs); -VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs); -VFPEvaluation operator*(double lhs, const VFPEvaluation& rhs); - -/** - * Helper function which interpolates data using the indices etc. given in the inputs. - */ -VFPEvaluation interpolate(const VFPProdTable& table, - const InterpData& flo_i, - const InterpData& thp_i, - const InterpData& wfr_i, - const InterpData& gfr_i, - const InterpData& alq_i); - -/** - * This basically models interpolate(VFPProdTable::array_type, ...) - * which performs 5D interpolation, but here for the 2D case only - */ -VFPEvaluation interpolate(const VFPInjTable& table, - const InterpData& flo_i, - const InterpData& thp_i); - -VFPEvaluation bhp(const VFPProdTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp, - const double alq, - const double explicit_wfr, - const double explicit_gfr, - const bool use_vfpexplicit); - -VFPEvaluation bhp(const VFPInjTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp); - - /** * Returns the table from the map if found, or throws an exception */ @@ -164,54 +127,95 @@ bool hasTable(const std::map>& tables, int return (entry != tables.end() ); } - /** * Returns the type variable for FLO/GFR/WFR for production tables */ template TYPE getType(const TABLE& table); +} -/** - * This function finds the value of THP given a specific BHP. - * Essentially: - * Given the function f(thp_array(x)) = bhp_array(x), which is piecewise linear, - * find thp so that f(thp) = bhp. - */ -double findTHP(const std::vector& bhp_array, - const std::vector& thp_array, - double bhp); +template +class VFPHelpers { +public: + /** + * Helper function to find indices etc. for linear interpolation and extrapolation + * @param value_in Value to find in values + * @param values Sorted list of values to search for value in. + * @return Data required to find the interpolated value + */ + static detail::InterpData findInterpData(const Scalar value_in, + const std::vector& values); -/** -* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq -*/ -std::pair -getMinimumBHPCoordinate(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq); + /** + * Helper function which interpolates data using the indices etc. given in the inputs. + */ + static detail::VFPEvaluation interpolate(const VFPProdTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i, + const detail::InterpData& wfr_i, + const detail::InterpData& gfr_i, + const detail::InterpData& alq_i); -/** -* Get (flo, bhp) at largest occuring stable vfp/ipr-intersection -* if it exists -*/ -std::optional> -intersectWithIPR(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq, - const double ipr_a, - const double ipr_b, - const std::function& adjust_bhp); + /** + * This basically models interpolate(VFPProdTable::array_type, ...) + * which performs 5D interpolation, but here for the 2D case only + */ + static detail::VFPEvaluation interpolate(const VFPInjTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i); -} // namespace detail + static detail::VFPEvaluation bhp(const VFPProdTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_vfpexplicit); + static detail::VFPEvaluation bhp(const VFPInjTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp); + + /** + * This function finds the value of THP given a specific BHP. + * Essentially: + * Given the function f(thp_array(x)) = bhp_array(x), which is piecewise linear, + * find thp so that f(thp) = bhp. + */ + static Scalar findTHP(const std::vector& bhp_array, + const std::vector& thp_array, + Scalar bhp); + + /** + * Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq + */ + static std::pair + getMinimumBHPCoordinate(const VFPProdTable& table, + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq); + + /** + * Get (flo, bhp) at largest occuring stable vfp/ipr-intersection + * if it exists + */ + static std::optional> + intersectWithIPR(const VFPProdTable& table, + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq, + const Scalar ipr_a, + const Scalar ipr_b, + const std::function& adjust_bhp); +}; } // namespace - - - #endif /* OPM_AUTODIFF_VFPHELPERS_HPP_ */ diff --git a/opm/simulators/wells/VFPInjProperties.cpp b/opm/simulators/wells/VFPInjProperties.cpp index 2da1e3749..0ba7be123 100644 --- a/opm/simulators/wells/VFPInjProperties.cpp +++ b/opm/simulators/wells/VFPInjProperties.cpp @@ -41,7 +41,7 @@ bhp(const int table_id, { const VFPInjTable& table = detail::getTable(m_tables, table_id); - detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg); + detail::VFPEvaluation retval = VFPHelpers::bhp(table, aqua, liquid, vapour, thp_arg); return retval.value; } @@ -61,7 +61,7 @@ thp(const int table_id, return 0.; } - const std::vector thp_array = table.getTHPAxis(); + const auto thp_array = table.getTHPAxis(); int nthp = thp_array.size(); /** @@ -69,14 +69,14 @@ thp(const int table_id, * by interpolating for every value of thp. This might be somewhat * expensive, but let us assome that nthp is small */ - auto flo_i = detail::findInterpData(flo, table.getFloAxis()); + const auto flo_i = VFPHelpers::findInterpData(flo, table.getFloAxis()); std::vector bhp_array(nthp); for (int i = 0; i < nthp; ++i) { - auto thp_i = detail::findInterpData(thp_array[i], thp_array); - bhp_array[i] = detail::interpolate(table, flo_i, thp_i).value; + auto thp_i = VFPHelpers::findInterpData(thp_array[i], thp_array); + bhp_array[i] = VFPHelpers::interpolate(table, flo_i, thp_i).value; } - return detail::findTHP(bhp_array, thp_array, bhp_arg); + return VFPHelpers::findTHP(bhp_array, thp_array, bhp_arg); } template @@ -115,10 +115,10 @@ EvalWell VFPInjProperties::bhp(const int table_id, //First, find the values to interpolate between //Value of FLO is negative in OPM for producers, but positive in VFP table - auto flo_i = detail::findInterpData(flo.value(), table.getFloAxis()); - auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant + const auto flo_i = VFPHelpers::findInterpData(flo.value(), table.getFloAxis()); + const auto thp_i = VFPHelpers::findInterpData(thp, table.getTHPAxis()); // assume constant - detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i); + detail::VFPEvaluation bhp_val = VFPHelpers::interpolate(table, flo_i, thp_i); bhp = bhp_val.dflo * flo; bhp.setValue(bhp_val.value); // thp is assumed constant i.e. diff --git a/opm/simulators/wells/VFPProdProperties.cpp b/opm/simulators/wells/VFPProdProperties.cpp index 8114d5b03..1ddec2d7c 100644 --- a/opm/simulators/wells/VFPProdProperties.cpp +++ b/opm/simulators/wells/VFPProdProperties.cpp @@ -67,17 +67,17 @@ thp(const int table_id, * by interpolating for every value of thp. This might be somewhat * expensive, but let us assome that nthp is small. */ - auto flo_i = detail::findInterpData( flo, table.getFloAxis()); - auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); - auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); - auto alq_i = detail::findInterpData( alq, table.getALQAxis()); + auto flo_i = VFPHelpers::findInterpData( flo, table.getFloAxis()); + auto wfr_i = VFPHelpers::findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = VFPHelpers::findInterpData( gfr, table.getGFRAxis()); + auto alq_i = VFPHelpers::findInterpData( alq, table.getALQAxis()); std::vector bhp_array(nthp); for (int i = 0; i < nthp; ++i) { - auto thp_i = detail::findInterpData(thp_array[i], thp_array); - bhp_array[i] = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; + auto thp_i = VFPHelpers::findInterpData(thp_array[i], thp_array); + bhp_array[i] = VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; } - return detail::findTHP(bhp_array, thp_array, bhp_arg); + return VFPHelpers::findTHP(bhp_array, thp_array, bhp_arg); } template @@ -94,7 +94,9 @@ bhp(const int table_id, { const VFPProdTable& table = detail::getTable(m_tables, table_id); - detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq, explicit_wfr,explicit_gfr, use_expvfp); + detail::VFPEvaluation retval = VFPHelpers::bhp(table, aqua, liquid, vapour, + thp_arg, alq, explicit_wfr, + explicit_gfr, use_expvfp); return retval.value; } @@ -124,16 +126,16 @@ bhpwithflo(const std::vector& flos, { // Get the table const VFPProdTable& table = detail::getTable(m_tables, table_id); - const auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant - const auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); - const auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); - const auto alq_i = detail::findInterpData( alq, table.getALQAxis()); //assume constant + const auto thp_i = VFPHelpers::findInterpData( thp, table.getTHPAxis()); // assume constant + const auto wfr_i = VFPHelpers::findInterpData( wfr, table.getWFRAxis()); + const auto gfr_i = VFPHelpers::findInterpData( gfr, table.getGFRAxis()); + const auto alq_i = VFPHelpers::findInterpData( alq, table.getALQAxis()); //assume constant std::vector bhps(flos.size(), 0.); for (std::size_t i = 0; i < flos.size(); ++i) { // Value of FLO is negative in OPM for producers, but positive in VFP table - const auto flo_i = detail::findInterpData(-flos[i], table.getFloAxis()); - const detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + const auto flo_i = VFPHelpers::findInterpData(-flos[i], table.getFloAxis()); + const detail::VFPEvaluation bhp_val = VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); // TODO: this kind of breaks the conventions for the functions here by putting dp within the function bhps[i] = bhp_val.value - dp; @@ -152,7 +154,7 @@ minimumBHP(const int table_id, { // Get the table const VFPProdTable& table = detail::getTable(m_tables, table_id); - const auto retval = detail::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq); + const auto retval = VFPHelpers::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq); // returned pair is (flo, bhp) return retval.second; } @@ -191,13 +193,14 @@ bhp(const int table_id, //First, find the values to interpolate between //Value of FLO is negative in OPM for producers, but positive in VFP table - auto flo_i = detail::findInterpData(-flo.value(), table.getFloAxis()); - auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant - auto wfr_i = detail::findInterpData( wfr.value(), table.getWFRAxis()); - auto gfr_i = detail::findInterpData( gfr.value(), table.getGFRAxis()); - auto alq_i = detail::findInterpData( alq, table.getALQAxis()); //assume constant + auto flo_i = VFPHelpers::findInterpData(-flo.value(), table.getFloAxis()); + auto thp_i = VFPHelpers::findInterpData( thp, table.getTHPAxis()); // assume constant + auto wfr_i = VFPHelpers::findInterpData( wfr.value(), table.getWFRAxis()); + auto gfr_i = VFPHelpers::findInterpData( gfr.value(), table.getGFRAxis()); + auto alq_i = VFPHelpers::findInterpData( alq, table.getALQAxis()); //assume constant - detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + detail::VFPEvaluation bhp_val = VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, + gfr_i, alq_i); bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo); bhp.setValue(bhp_val.value); diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 14a4cabaf..b90ddf8e9 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -919,7 +919,9 @@ isStableSolution(const WellState& well_state, const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); const bool use_vfpexplicit = well_.useVfpExplicit(); - detail::VFPEvaluation bhp = detail::bhp(table, aqua, liquid, vapour, thp, well_.getALQ(well_state), wfr, gfr, use_vfpexplicit); + auto bhp = VFPHelpers::bhp(table, aqua, liquid, vapour, thp, + well_.getALQ(well_state), wfr, gfr, + use_vfpexplicit); if (bhp.dflo >= 0) { return true; @@ -964,7 +966,10 @@ estimateStableBhp(const WellState& well_state, auto bhp_adjusted = [this, &thp, &dp_hydro](const Scalar bhp) { return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp); }; - const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted); + const auto retval = VFPHelpers::intersectWithIPR(table, thp, wfr, gfr, + well_.getALQ(well_state), + ipr.first, ipr.second, + bhp_adjusted); if (retval.has_value()) { // returned pair is (flo, bhp) return retval.value().second; diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 44e17b047..fac8740a7 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -68,12 +68,12 @@ BOOST_AUTO_TEST_CASE(findInterpData) double first = 1; double last = 15; - Opm::detail::InterpData eval0 = Opm::detail::findInterpData(exact, values); - Opm::detail::InterpData eval1 = Opm::detail::findInterpData(interpolate, values); - Opm::detail::InterpData eval2 = Opm::detail::findInterpData(extrapolate_left, values); - Opm::detail::InterpData eval3 = Opm::detail::findInterpData(extrapolate_right, values); - Opm::detail::InterpData eval4 = Opm::detail::findInterpData(first, values); - Opm::detail::InterpData eval5 = Opm::detail::findInterpData(last, values); + auto eval0 = Opm::VFPHelpers::findInterpData(exact, values); + auto eval1 = Opm::VFPHelpers::findInterpData(interpolate, values); + auto eval2 = Opm::VFPHelpers::findInterpData(extrapolate_left, values); + auto eval3 = Opm::VFPHelpers::findInterpData(extrapolate_right, values); + auto eval4 = Opm::VFPHelpers::findInterpData(first, values); + auto eval5 = Opm::VFPHelpers::findInterpData(last, values); BOOST_CHECK_EQUAL(eval0.ind_[0], 2); BOOST_CHECK_EQUAL(eval0.ind_[1], 3); @@ -109,7 +109,7 @@ BOOST_AUTO_TEST_SUITE_END() // HelperTests * values data is given at */ struct TrivialFixture { - typedef Opm::detail::VFPEvaluation VFPEvaluation; + typedef Opm::detail::VFPEvaluation VFPEvaluation; TrivialFixture() : table_ids(1, 1), thp_axis{0.0, 1.0},