OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-07-07 04:53:03 -05:00
Code Issues Packages Projects Releases Wiki Activity

Use std::reference_wrapper for VFP tables

Browse Source
This commit is contained in:
Joakim Hove 2021-01-25 10:08:31 +01:00
parent 57d158bbbe
commit e4789d4eb7
11 changed files with 117 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/printvfp.cpp
View File

@ -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) {

4
opm/simulators/wells/GasLiftRuntime_impl.hpp
View File

@ -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:

12
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -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);

12
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -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);

54
opm/simulators/wells/VFPHelpers.hpp
View File

@ -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 <typename T>
const T* getTable(const std::map<int, T*> tables, int table_id) {
const T& getTable(const std::map<int, std::reference_wrapper<const T>> 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<int, T*> tables, int table_id) {
* Check whether we have a table with the table number
*/
template <typename T>
bool hasTable(const std::map<int, T*> tables, int table_id) {
bool hasTable(const std::map<int, std::reference_wrapper<const T>> tables, int table_id) {
const auto entry = tables.find(table_id);
return (entry != tables.end() );
}
@ -584,24 +584,24 @@ bool hasTable(const std::map<int, T*> tables, int table_id) {
* Returns the type variable for FLO/GFR/WFR for production tables
*/
template <typename TYPE, typename TABLE>
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();
}

18
opm/simulators/wells/VFPInjProperties.cpp
View File

@ -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<double> thp_array = table->getTHPAxis();
const std::vector<double> 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<double> 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;
bhp_array[i] = detail::interpolate(table, flo_i, thp_i).value;
}
double retval = detail::findTHP(bhp_array, thp_array, bhp_arg);
return retval;
}
const VFPInjTable* VFPInjProperties::getTable(const int table_id) const {
const VFPInjTable& VFPInjProperties::getTable(const int table_id) const {
return detail::getTable(m_tables, table_id);
}
@ -83,8 +83,8 @@ bool VFPInjProperties::hasTable(const int table_id) const {
return detail::hasTable(m_tables, table_id);
}
void VFPInjProperties::addTable(const VFPInjTable * new_table) {
this->m_tables.emplace( new_table->getTableNum(), new_table );
void VFPInjProperties::addTable(const VFPInjTable& new_table) {
this->m_tables.emplace( new_table.getTableNum(), new_table );
}

30
opm/simulators/wells/VFPInjProperties.hpp
View File

@ -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<int, const VFPInjTable*> m_tables;
std::map<int, std::reference_wrapper<const VFPInjTable>> m_tables;
};

56
opm/simulators/wells/VFPProdProperties.cpp
View File

@ -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<double> thp_array = table->getTHPAxis();
const std::vector<double> 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<double> 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;
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);
@ -87,14 +87,14 @@ double VFPProdProperties::bhp(int table_id,
const double& vapour,
const double& thp_arg,
const double& alq) const {
const VFPProdTable* table = detail::getTable(m_tables, table_id);
const VFPProdTable& table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq);
return retval.value;
}
const VFPProdTable* VFPProdProperties::getTable(const int table_id) const {
const VFPProdTable& VFPProdProperties::getTable(const int table_id) const {
return detail::getTable(m_tables, table_id);
}
@ -114,17 +114,17 @@ bhpwithflo(const std::vector<double>& 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<double> 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<double>& 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<double> flo_samples = table->getFloAxis();
std::vector<double> flo_samples = table.getFloAxis();
if (flo_samples[0] > 0.) {
flo_samples.insert(flo_samples.begin(), 0.);
@ -244,8 +244,8 @@ calculateBhpWithTHPTarget(const std::vector<double>& 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 );
}
}

40
opm/simulators/wells/VFPProdProperties.hpp
View File

@ -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<int, const VFPProdTable*> m_tables;
std::map<int, std::reference_wrapper<const VFPProdTable>> m_tables;
};

12
opm/simulators/wells/VFPProperties.hpp
View File

@ -46,14 +46,14 @@ public:
* @param prod_tables A map of different VFPPROD tables.
*/
VFPProperties(const std::vector<const VFPInjTable *>& inj_tables,
const std::vector<const VFPProdTable *>& prod_tables)
VFPProperties(const std::vector<std::reference_wrapper<const VFPInjTable>>& inj_tables,
const std::vector<std::reference_wrapper<const VFPProdTable>>& 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 );
};
/**

6
tests/test_vfpproperties.cpp
View File

@ -206,7 +206,7 @@ struct TrivialFixture {
data));
properties = std::make_shared<Opm::VFPProdProperties>();
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