From e4789d4eb7617332aa7532978e805026c93b2757 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Mon, 25 Jan 2021 10:08:31 +0100 Subject: [PATCH] Use std::reference_wrapper for VFP tables --- examples/printvfp.cpp | 2 +- opm/simulators/wells/GasLiftRuntime_impl.hpp | 4 +- .../wells/MultisegmentWell_impl.hpp | 12 ++-- opm/simulators/wells/StandardWell_impl.hpp | 12 ++-- opm/simulators/wells/VFPHelpers.hpp | 54 +++++++++--------- opm/simulators/wells/VFPInjProperties.cpp | 18 +++--- opm/simulators/wells/VFPInjProperties.hpp | 30 ++++------ opm/simulators/wells/VFPProdProperties.cpp | 56 +++++++++---------- opm/simulators/wells/VFPProdProperties.hpp | 40 ++++++------- opm/simulators/wells/VFPProperties.hpp | 12 ++-- tests/test_vfpproperties.cpp | 6 +- 11 files changed, 117 insertions(+), 129 deletions(-) diff --git a/examples/printvfp.cpp b/examples/printvfp.cpp index b29e81de1..b52199948 100644 --- a/examples/printvfp.cpp +++ b/examples/printvfp.cpp @@ -115,7 +115,7 @@ int main(int argc, char** argv) thps[ii] = (1.0 - q) * thp_min + q * thp_max; } - const VFPProdTable& table = *(setup.vfp_properties->getProd()->getTable(table_id)); + const VFPProdTable& table = setup.vfp_properties->getProd()->getTable(table_id); std::cout.precision(12); for (double rate : rates) { for (double thp : thps) { diff --git a/opm/simulators/wells/GasLiftRuntime_impl.hpp b/opm/simulators/wells/GasLiftRuntime_impl.hpp index e3115dcae..165e2a40f 100644 --- a/opm/simulators/wells/GasLiftRuntime_impl.hpp +++ b/opm/simulators/wells/GasLiftRuntime_impl.hpp @@ -420,8 +420,8 @@ setAlqMaxRate_(const GasLiftOpt::Well &well) // According to the Eclipse manual for WLIFTOPT, item 3: // The default value should be set to the largest ALQ // value in the well's VFP table - const auto& table = *(std_well_.vfp_properties_->getProd()->getTable( - this->controls_.vfp_table_number)); + const auto& table = std_well_.vfp_properties_->getProd()->getTable( + this->controls_.vfp_table_number); const auto& alq_values = table.getALQAxis(); // Assume the alq_values are sorted in ascending order, so // the last item should be the largest value: diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index d55d4f54d..0d243dfde 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -2144,7 +2144,7 @@ namespace Opm double thp = 0.0; if (this->isInjector()) { const int table_id = well_ecl_.vfp_table_number(); - const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); @@ -2152,7 +2152,7 @@ namespace Opm else if (this->isProducer()) { const int table_id = well_ecl_.vfp_table_number(); const double alq = well_ecl_.alq_value(); - const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); thp = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq); @@ -2193,13 +2193,13 @@ namespace Opm if (well.isInjector() ) { const auto& controls = well.injectionControls(summaryState); - const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState)) - dp; } else if (well.isProducer()) { const auto& controls = well.productionControls(summaryState); - const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState), controls.alq_value) - dp; } @@ -3605,7 +3605,7 @@ namespace Opm // Make the fbhp() function. const auto& controls = well_ecl_.productionControls(summary_state); - const auto& table = *(vfp_properties_->getProd()->getTable(controls.vfp_table_number)); + const auto& table = vfp_properties_->getProd()->getTable(controls.vfp_table_number); const double vfp_ref_depth = table.getDatumDepth(); const double rho = segment_densities_[0].value(); // Use the density at the top perforation. const double thp_limit = this->getTHPConstraint(summary_state); @@ -3829,7 +3829,7 @@ namespace Opm // Make the fbhp() function. const auto& controls = well_ecl_.injectionControls(summary_state); - const auto& table = *(vfp_properties_->getInj()->getTable(controls.vfp_table_number)); + const auto& table = vfp_properties_->getInj()->getTable(controls.vfp_table_number); const double vfp_ref_depth = table.getDatumDepth(); const double rho = segment_densities_[0].value(); // Use the density at the top perforation. const double thp_limit = this->getTHPConstraint(summary_state); diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 9906fff8b..026f69fd3 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -3027,13 +3027,13 @@ namespace Opm if (this->isInjector() ) { const auto& controls = well.injectionControls(summaryState); - const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState)) - dp; } else if (this->isProducer()) { const auto& controls = well.productionControls(summaryState); - const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState), getALQ(well_state)) - dp; } @@ -3066,7 +3066,7 @@ namespace Opm double thp = 0.0; if (this->isInjector()) { const int table_id = well_ecl_.vfp_table_number(); - const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); @@ -3074,7 +3074,7 @@ namespace Opm else if (this->isProducer()) { const int table_id = well_ecl_.vfp_table_number(); const double alq = getALQ(well_state); - const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id).getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); thp = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq); @@ -3714,7 +3714,7 @@ namespace Opm // Make the fbhp() function. const auto& controls = well_ecl_.productionControls(summary_state); - const auto& table = *(vfp_properties_->getProd()->getTable(controls.vfp_table_number)); + const auto& table = vfp_properties_->getProd()->getTable(controls.vfp_table_number); const double vfp_ref_depth = table.getDatumDepth(); const double rho = perf_densities_[0]; // Use the density at the top perforation. const double thp_limit = this->getTHPConstraint(summary_state); @@ -3920,7 +3920,7 @@ namespace Opm // Make the fbhp() function. const auto& controls = well_ecl_.injectionControls(summary_state); - const auto& table = *(vfp_properties_->getInj()->getTable(controls.vfp_table_number)); + const auto& table = vfp_properties_->getInj()->getTable(controls.vfp_table_number); const double vfp_ref_depth = table.getDatumDepth(); const double rho = perf_densities_[0]; // Use the density at the top perforation. const double thp_limit = this->getTHPConstraint(summary_state); diff --git a/opm/simulators/wells/VFPHelpers.hpp b/opm/simulators/wells/VFPHelpers.hpp index 60a965df6..87a0a6f5a 100644 --- a/opm/simulators/wells/VFPHelpers.hpp +++ b/opm/simulators/wells/VFPHelpers.hpp @@ -503,26 +503,26 @@ inline VFPEvaluation interpolate( return nn[0][0]; } -inline VFPEvaluation bhp(const VFPProdTable* table, +inline VFPEvaluation bhp(const VFPProdTable& table, const double& aqua, const double& liquid, const double& vapour, const double& thp, const double& alq) { //Find interpolation variables - double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); - double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType()); - double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); + double flo = detail::getFlo(aqua, liquid, vapour, table.getFloType()); + double wfr = detail::getWFR(aqua, liquid, vapour, table.getWFRType()); + double gfr = detail::getGFR(aqua, liquid, vapour, table.getGFRType()); //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 = 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()); - detail::VFPEvaluation retval = detail::interpolate(*table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); return retval; } @@ -531,20 +531,20 @@ inline VFPEvaluation bhp(const VFPProdTable* table, -inline VFPEvaluation bhp(const VFPInjTable* table, +inline VFPEvaluation bhp(const VFPInjTable& table, const double& aqua, const double& liquid, const double& vapour, const double& thp) { //Find interpolation variables - double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + double flo = detail::getFlo(aqua, liquid, vapour, table.getFloType()); //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 = detail::findInterpData(flo, table.getFloAxis()); + auto thp_i = detail::findInterpData(thp, table.getTHPAxis()); //Then perform the interpolation itself - detail::VFPEvaluation retval = detail::interpolate(*table, flo_i, thp_i); + detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i); return retval; } @@ -560,13 +560,13 @@ inline VFPEvaluation bhp(const VFPInjTable* table, * Returns the table from the map if found, or throws an exception */ template -const T* getTable(const std::map tables, int table_id) { +const T& getTable(const std::map> tables, int table_id) { auto entry = tables.find(table_id); if (entry == tables.end()) { OPM_THROW(std::invalid_argument, "Nonexistent VFP table " << table_id << " referenced."); } else { - return entry->second; + return entry->second.get(); } } @@ -574,7 +574,7 @@ const T* getTable(const std::map tables, int table_id) { * Check whether we have a table with the table number */ template -bool hasTable(const std::map tables, int table_id) { +bool hasTable(const std::map> tables, int table_id) { const auto entry = tables.find(table_id); return (entry != tables.end() ); } @@ -584,24 +584,24 @@ bool hasTable(const std::map tables, int table_id) { * Returns the type variable for FLO/GFR/WFR for production tables */ template -TYPE getType(const TABLE* table); +TYPE getType(const TABLE& table); template <> inline -VFPProdTable::FLO_TYPE getType(const VFPProdTable* table) { - return table->getFloType(); +VFPProdTable::FLO_TYPE getType(const VFPProdTable& table) { + return table.getFloType(); } template <> inline -VFPProdTable::WFR_TYPE getType(const VFPProdTable* table) { - return table->getWFRType(); +VFPProdTable::WFR_TYPE getType(const VFPProdTable& table) { + return table.getWFRType(); } template <> inline -VFPProdTable::GFR_TYPE getType(const VFPProdTable* table) { - return table->getGFRType(); +VFPProdTable::GFR_TYPE getType(const VFPProdTable& table) { + return table.getGFRType(); } @@ -610,8 +610,8 @@ VFPProdTable::GFR_TYPE getType(const VFPProdTable* table) { */ template <> inline -VFPInjTable::FLO_TYPE getType(const VFPInjTable* table) { - return table->getFloType(); +VFPInjTable::FLO_TYPE getType(const VFPInjTable& table) { + return table.getFloType(); } diff --git a/opm/simulators/wells/VFPInjProperties.cpp b/opm/simulators/wells/VFPInjProperties.cpp index 90be1980b..6a2863acb 100644 --- a/opm/simulators/wells/VFPInjProperties.cpp +++ b/opm/simulators/wells/VFPInjProperties.cpp @@ -39,7 +39,7 @@ double VFPInjProperties::bhp(int table_id, const double& liquid, const double& vapour, const double& thp_arg) const { - const VFPInjTable* table = detail::getTable(m_tables, table_id); + const VFPInjTable& table = detail::getTable(m_tables, table_id); detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg); return retval.value; @@ -51,12 +51,12 @@ double VFPInjProperties::thp(int table_id, const double& liquid, const double& vapour, const double& bhp_arg) const { - const VFPInjTable* table = detail::getTable(m_tables, table_id); + const VFPInjTable& table = detail::getTable(m_tables, table_id); //Find interpolation variables - double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + double flo = detail::getFlo(aqua, liquid, vapour, table.getFloType()); - const std::vector thp_array = table->getTHPAxis(); + const std::vector thp_array = table.getTHPAxis(); int nthp = thp_array.size(); /** @@ -64,18 +64,18 @@ double VFPInjProperties::thp(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 flo_i = detail::findInterpData(flo, table.getFloAxis()); std::vector bhp_array(nthp); for (int i=0; im_tables.emplace( new_table->getTableNum(), new_table ); +void VFPInjProperties::addTable(const VFPInjTable& new_table) { + this->m_tables.emplace( new_table.getTableNum(), new_table ); } diff --git a/opm/simulators/wells/VFPInjProperties.hpp b/opm/simulators/wells/VFPInjProperties.hpp index 9f985ac6f..8a14010fc 100644 --- a/opm/simulators/wells/VFPInjProperties.hpp +++ b/opm/simulators/wells/VFPInjProperties.hpp @@ -38,7 +38,7 @@ public: /** * Takes *no* ownership of data. */ - void addTable(const VFPInjTable * new_table); + void addTable(const VFPInjTable& new_table); /** * Linear interpolation of bhp as a function of the input parameters given as @@ -64,27 +64,21 @@ public: const double& thp) const { //Get the table - const VFPInjTable* table = detail::getTable(m_tables, table_id); + const VFPInjTable& table = detail::getTable(m_tables, table_id); EvalWell bhp = 0.0 * aqua; //Find interpolation variables - EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + EvalWell flo = detail::getFlo(aqua, liquid, vapour, table.getFloType()); - //Compute the BHP for each well independently - if (table != nullptr) { - //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 + //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 - detail::VFPEvaluation bhp_val = detail::interpolate(*table, flo_i, thp_i); + detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i); - bhp = bhp_val.dflo * flo; - bhp.setValue(bhp_val.value); // thp is assumed constant i.e. - } - else { - bhp.setValue(-1e100); //Signal that this value has not been calculated properly, due to "missing" table - } + bhp = bhp_val.dflo * flo; + bhp.setValue(bhp_val.value); // thp is assumed constant i.e. return bhp; } @@ -92,7 +86,7 @@ public: * Returns the table associated with the ID, or throws an exception if * the table does not exist */ - const VFPInjTable* getTable(const int table_id) const; + const VFPInjTable& getTable(const int table_id) const; /** * Check whether there is table associated with ID @@ -142,7 +136,7 @@ public: protected: // Map which connects the table number with the table itself - std::map m_tables; + std::map> m_tables; }; diff --git a/opm/simulators/wells/VFPProdProperties.cpp b/opm/simulators/wells/VFPProdProperties.cpp index 168c3c35f..41d369f06 100644 --- a/opm/simulators/wells/VFPProdProperties.cpp +++ b/opm/simulators/wells/VFPProdProperties.cpp @@ -39,7 +39,7 @@ double VFPProdProperties::thp(int table_id, const double& vapour, const double& bhp_arg, const double& alq) const { - const VFPProdTable* table = detail::getTable(m_tables, table_id); + const VFPProdTable& table = detail::getTable(m_tables, table_id); // Find interpolation variables. double flo = 0.0; @@ -49,16 +49,16 @@ double VFPProdProperties::thp(int table_id, // All zero, likely at initial state. // Set FLO variable to minimum to avoid extrapolation. // The water and gas fractions are kept at 0.0. - flo = table->getFloAxis().front(); + flo = table.getFloAxis().front(); } else { // The usual case. // Recall that production rate is negative in Opm, so switch the sign. - flo = -detail::getFlo(aqua, liquid, vapour, table->getFloType()); - wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType()); - gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); + flo = -detail::getFlo(aqua, liquid, vapour, table.getFloType()); + wfr = detail::getWFR(aqua, liquid, vapour, table.getWFRType()); + gfr = detail::getGFR(aqua, liquid, vapour, table.getGFRType()); } - const std::vector thp_array = table->getTHPAxis(); + const std::vector thp_array = table.getTHPAxis(); int nthp = thp_array.size(); /** @@ -66,14 +66,14 @@ double VFPProdProperties::thp(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 = 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()); std::vector bhp_array(nthp); for (int i=0; i& flos, const double dp) const { // 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 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 std::vector bhps(flos.size(), 0.); for (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 = detail::findInterpData(-flos[i], table.getFloAxis()); + const detail::VFPEvaluation bhp_val = detail::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; @@ -177,24 +177,24 @@ calculateBhpWithTHPTarget(const std::vector& ipr_a, const int Oil = BlackoilPhases::Liquid; const int Gas = BlackoilPhases::Vapour; - const VFPProdTable* table = detail::getTable(m_tables, thp_table_id); + const VFPProdTable& table = detail::getTable(m_tables, thp_table_id); const double aqua_bhp_limit = rates_bhp_limit[Water]; const double liquid_bhp_limit = rates_bhp_limit[Oil]; const double vapour_bhp_limit = rates_bhp_limit[Gas]; - const double flo_bhp_limit = detail::getFlo(aqua_bhp_limit, liquid_bhp_limit, vapour_bhp_limit, table->getFloType() ); + const double flo_bhp_limit = detail::getFlo(aqua_bhp_limit, liquid_bhp_limit, vapour_bhp_limit, table.getFloType() ); const double aqua_bhp_middle = rates_bhp_middle[Water]; const double liquid_bhp_middle = rates_bhp_middle[Oil]; const double vapour_bhp_middle = rates_bhp_middle[Gas]; - const double flo_bhp_middle = detail::getFlo(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table->getFloType() ); + const double flo_bhp_middle = detail::getFlo(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table.getFloType() ); // we use the ratios based on the middle value of bhp_limit and bhp_safe_limit - const double wfr = detail::getWFR(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table->getWFRType()); - const double gfr = detail::getGFR(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table->getGFRType()); + const double wfr = detail::getWFR(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table.getWFRType()); + const double gfr = detail::getGFR(aqua_bhp_middle, liquid_bhp_middle, vapour_bhp_middle, table.getGFRType()); // we get the flo sampling points from the table, // then extend it with zero and rate under bhp_limit for extrapolation - std::vector flo_samples = table->getFloAxis(); + std::vector flo_samples = table.getFloAxis(); if (flo_samples[0] > 0.) { flo_samples.insert(flo_samples.begin(), 0.); @@ -244,8 +244,8 @@ calculateBhpWithTHPTarget(const std::vector& ipr_a, } } -void VFPProdProperties::addTable(const VFPProdTable * new_table) { - this->m_tables.emplace( new_table->getTableNum(), new_table ); +void VFPProdProperties::addTable(const VFPProdTable& new_table) { + this->m_tables.emplace( new_table.getTableNum(), new_table ); } } diff --git a/opm/simulators/wells/VFPProdProperties.hpp b/opm/simulators/wells/VFPProdProperties.hpp index 6a60486f3..153fd0d5f 100644 --- a/opm/simulators/wells/VFPProdProperties.hpp +++ b/opm/simulators/wells/VFPProdProperties.hpp @@ -42,7 +42,7 @@ public: /** * Takes *no* ownership of data. */ - void addTable(const VFPProdTable * new_table); + void addTable(const VFPProdTable& new_table); /** * Linear interpolation of bhp as a function of the input parameters given as @@ -70,32 +70,26 @@ public: const double& alq) const { //Get the table - const VFPProdTable* table = detail::getTable(m_tables, table_id); + const VFPProdTable& table = detail::getTable(m_tables, table_id); EvalWell bhp = 0.0 * aqua; //Find interpolation variables - EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); - EvalWell wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType()); - EvalWell gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); + EvalWell flo = detail::getFlo(aqua, liquid, vapour, table.getFloType()); + EvalWell wfr = detail::getWFR(aqua, liquid, vapour, table.getWFRType()); + EvalWell gfr = detail::getGFR(aqua, liquid, vapour, table.getGFRType()); - //Compute the BHP for each well independently - if (table != nullptr) { - //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 + //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 - detail::VFPEvaluation bhp_val = detail::interpolate(*table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); - bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo * flo); - bhp.setValue(bhp_val.value); - } - else { - bhp.setValue(-1e100); //Signal that this value has not been calculated properly, due to "missing" table - } + bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo * flo); + bhp.setValue(bhp_val.value); return bhp; } @@ -141,7 +135,7 @@ public: * Returns the table associated with the ID, or throws an exception if * the table does not exist */ - const VFPProdTable* getTable(const int table_id) const; + const VFPProdTable& getTable(const int table_id) const; /** * Check whether there is table associated with ID @@ -180,7 +174,7 @@ protected: const double dp) const; // Map which connects the table number with the table itself - std::map m_tables; + std::map> m_tables; }; diff --git a/opm/simulators/wells/VFPProperties.hpp b/opm/simulators/wells/VFPProperties.hpp index 96c5661c1..00586a92e 100644 --- a/opm/simulators/wells/VFPProperties.hpp +++ b/opm/simulators/wells/VFPProperties.hpp @@ -46,14 +46,14 @@ public: * @param prod_tables A map of different VFPPROD tables. */ - VFPProperties(const std::vector& inj_tables, - const std::vector& prod_tables) + VFPProperties(const std::vector>& inj_tables, + const std::vector>& prod_tables) { - for (const auto& vfpinj_ptr : inj_tables) - this->m_inj.addTable( vfpinj_ptr ); + for (const auto& vfpinj : inj_tables) + this->m_inj.addTable( vfpinj ); - for (const auto& vfpprod_ptr : prod_tables) - this->m_prod.addTable( vfpprod_ptr ); + for (const auto& vfpprod : prod_tables) + this->m_prod.addTable( vfpprod ); }; /** diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 11a133659..b62e4944f 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -206,7 +206,7 @@ struct TrivialFixture { data)); properties = std::make_shared(); - properties->addTable( table.get() ); + properties->addTable( *table ); } double& operator()(size_t thp_idx, size_t wfr_idx, size_t gfr_idx, size_t alq_idx, size_t flo_idx) { @@ -592,7 +592,7 @@ VFPPROD \n\ auto deck = parser.parseString(table_str); Opm::VFPProdTable table(deck.getKeyword("VFPPROD", 0), units); Opm::VFPProdProperties properties; - properties.addTable( &table ); + properties.addTable( table ); const int n = 5; //Number of points to check per axis double bhp_sad = 0.0; //Sum of absolute difference @@ -652,7 +652,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD) Opm::VFPProdTable table(deck.getKeyword("VFPPROD", 0), units); Opm::VFPProdProperties properties; - properties.addTable(&table); + properties.addTable(table); //Do some rudimentary testing //Get the BHP as a function of rate, thp, wfr, gfr, alq