diff --git a/opm/parser/eclipse/EclipseState/EclipseState.cpp b/opm/parser/eclipse/EclipseState/EclipseState.cpp index 982cd5c15..29383c7a5 100644 --- a/opm/parser/eclipse/EclipseState/EclipseState.cpp +++ b/opm/parser/eclipse/EclipseState/EclipseState.cpp @@ -246,6 +246,10 @@ namespace Opm { return m_watvisctTables; } + const std::vector& EclipseState::getVFPProdTables() const { + return m_vfpprodTables; + } + const std::vector& EclipseState::getSgofTables() const { return m_sgofTables; } @@ -357,6 +361,8 @@ namespace Opm { initPlyshlogTables(deck, "PLYSHLOG", m_plyshlogTables); + initVFPProdTables(deck, "VFPPROD", m_vfpprodTables); + // the ROCKTAB table comes with additional fun because the number of columns //depends on the presence of the RKTRMDIR keyword... initRocktabTables(deck); @@ -652,6 +658,25 @@ namespace Opm { } + void EclipseState::initVFPProdTables(DeckConstPtr deck, + const std::string& keywordName, + std::vector& tableVector) { + if (!deck->hasKeyword(keywordName)) { + return; + } + + if (!deck->numKeywords(keywordName)) { + complainAboutAmbiguousKeyword(deck, keywordName); + return; + } + + const auto& keyword = deck->getKeyword(keywordName); + + tableVector.push_back(VFPProdTable()); + + tableVector[0].init(keyword); + } + bool EclipseState::supportsGridProperty(const std::string& keyword, int enabledTypes) const { bool result = false; if (enabledTypes & IntProperties) diff --git a/opm/parser/eclipse/EclipseState/EclipseState.hpp b/opm/parser/eclipse/EclipseState/EclipseState.hpp index 55e6d900f..b6cdfc2aa 100644 --- a/opm/parser/eclipse/EclipseState/EclipseState.hpp +++ b/opm/parser/eclipse/EclipseState/EclipseState.hpp @@ -66,6 +66,7 @@ #include #include #include +#include "Tables/VFPProdTable.hpp" namespace Opm { class EclipseState { @@ -135,6 +136,7 @@ namespace Opm { const std::vector& getSwofTables() const; const std::vector& getSwfnTables() const; const std::vector& getWatvisctTables() const; + const std::vector& getVFPProdTables() const; size_t getNumPhases() const; // the unit system used by the deck. note that it is rarely needed to convert @@ -220,6 +222,9 @@ namespace Opm { void initPlyshlogTables(DeckConstPtr deck, const std::string& keywordName, std::vector& tableVector); + void initVFPProdTables(DeckConstPtr deck, + const std::string& keywordName, + std::vector& tableVector); void setMULTFLT(std::shared_ptr section) const; @@ -281,6 +286,7 @@ namespace Opm { std::vector m_swofTables; std::vector m_swfnTables; std::vector m_watvisctTables; + std::vector m_vfpprodTables; std::set phases; std::string m_title; diff --git a/opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp b/opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp new file mode 100644 index 000000000..3278f17b7 --- /dev/null +++ b/opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp @@ -0,0 +1,536 @@ +/* + Copyright 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_PARSER_ECLIPSE_ECLIPSESTATE_TABLES_VFPPRODTABLE_HPP_ +#define OPM_PARSER_ECLIPSE_ECLIPSESTATE_TABLES_VFPPRODTABLE_HPP_ + +#include +#include + + + +namespace Opm { + +/** + * Class for reading data from a VFPPROD (vertical flow performance production) table + */ +class VFPProdTable { +public: + typedef boost::multi_array array_type; + typedef boost::array extents; + + ///Rate type + enum FLO_TYPE { + FLO_OIL=1, //< Oil rate + FLO_LIQ, //< Liquid rate + FLO_GAS, //< Gas rate + FLO_INVALID + }; + + ///Water fraction variable + enum WFR_TYPE { + WFR_WOR=11, //< Water-oil ratio + WFR_WCT, //< Water cut + WFR_WGR, //< Water-gas ratio + WFR_INVALID + }; + + ///Gas fraction variable + enum GFR_TYPE { + GFR_GOR=21, //< Gas-oil ratio + GFR_GLR, //< Gas-liquid ratio + GFR_OGR, //< Oil-gas ratio + GFR_INVALID + }; + + ///Artificial lift quantity + enum ALQ_TYPE { + ALQ_GRAT=31, //< Lift as injection rate + ALQ_IGLR, //< Injection gas-liquid ratio + ALQ_TGLR, //< Total gas-liquid ratio + ALQ_PUMP, //< Pump rating + ALQ_COMP, //< Compressor power + ALQ_BEAN, //< Choke diameter + ALQ_UNDEF, //< Undefined + ALQ_INVALID + }; + + /** + * Constructor + */ + VFPProdTable() : m_table_num(-1), + m_datum_depth(-1), + m_flo_type(FLO_INVALID), + m_wfr_type(WFR_INVALID), + m_gfr_type(GFR_INVALID), + m_alq_type(ALQ_INVALID) { + + } + + /** + * Initializes objects from raw data. NOTE: All raw data assumed to be in SI units + * @param table_num VFP table number + * @param datum_depth Reference depth for BHP + * @param flo_type Specifies what flo_data represents + * @param wfr_type Specifies what wfr_data represents + * @param gfr_type Specifies what gfr_data represents + * @param alq_type Specifies what alq_data represents + * @param flo_data Axis for flo_type + * @param thp_data Axis for thp_type + * @param wfr_data Axis for wfr_type + * @param gfr_data Axis for gfr_type + * @param alq_data Axis for alq_type + * @param data BHP to be interpolated. Given as a 5D array so that + * BHP = data[thp][wfr][gfr][alq][flo] for the indices thp, wfr, etc. + */ + void init(int table_num, + double datum_depth, + FLO_TYPE flo_type, + WFR_TYPE wfr_type, + GFR_TYPE gfr_type, + ALQ_TYPE alq_type, + const std::vector& flo_data, + const std::vector& thp_data, + const std::vector& wfr_data, + const std::vector& gfr_data, + const std::vector& alq_data, + array_type data) { + m_table_num = table_num; + m_datum_depth = datum_depth; + m_flo_type = flo_type; + m_wfr_type = wfr_type; + m_gfr_type = gfr_type; + m_alq_type = alq_type; + m_flo_data = flo_data; + m_thp_data = thp_data; + m_wfr_data = wfr_data; + m_gfr_data = gfr_data; + m_alq_data = alq_data; + m_data = data; + check(); + } + + /** + * Constructor which parses a deck keyword and retrieves the relevant parts for a + * VFP table. + */ + void init(DeckKeywordConstPtr table) { + auto iter = table->begin(); + + auto header = (*iter++); + + assert(itemValid(header, "TABLE")); + m_table_num = header->getItem("TABLE")->getInt(0); + + assert(itemValid(header, "DATUM_DEPTH")); + m_datum_depth = header->getItem("DATUM_DEPTH")->getRawDouble(0); + + //Rate type + assert(itemValid(header, "RATE_TYPE")); + std::string flo_string = header->getItem("RATE_TYPE")->getString(0); + if (flo_string == "OIL") { + m_flo_type = FLO_OIL; + } + else if (flo_string == "LIQ") { + m_flo_type = FLO_LIQ; + } + else if (flo_string == "GAS") { + m_flo_type = FLO_GAS; + } + else { + m_flo_type = FLO_INVALID; + throw std::invalid_argument("Invalid RATE_TYPE string"); + } + + //Water fraction + assert(itemValid(header, "WFR")); + std::string wfr_string = header->getItem("WFR")->getString(0); + if (wfr_string == "WOR") { + m_wfr_type = WFR_WOR; + } + else if (wfr_string == "WCT") { + m_wfr_type = WFR_WCT; + } + else if (wfr_string == "WGR") { + m_wfr_type = WFR_WGR; + } + else { + m_wfr_type = WFR_INVALID; + throw std::invalid_argument("Invalid WFR string"); + } + + //Gas fraction + assert(itemValid(header, "GFR")); + std::string gfr_string = header->getItem("GFR")->getString(0); + if (gfr_string == "GOR") { + m_gfr_type = GFR_GOR; + } + else if (gfr_string == "GLR") { + m_gfr_type = GFR_GLR; + } + else if (gfr_string == "OGR") { + m_gfr_type = GFR_OGR; + } + else { + m_gfr_type = GFR_INVALID; + throw std::invalid_argument("Invalid GFR string"); + } + + //Definition of THP values, must be THP + if (itemValid(header, "PRESSURE_DEF")) { + std::string quantity_string = header->getItem("PRESSURE_DEF")->getString(0); + assert(quantity_string == "THP"); + } + + //Artificial lift + if (itemValid(header, "ALQ_DEF")) { + std::string alq_string = header->getItem("ALQ_DEF")->getString(0); + if (alq_string == "GRAT") { + m_alq_type = ALQ_GRAT; + } + else if (alq_string == "IGLR") { + m_alq_type = ALQ_IGLR; + } + else if (alq_string == "TGLR") { + m_alq_type = ALQ_TGLR; + } + else if (alq_string == "PUMP") { + m_alq_type = ALQ_PUMP; + } + else if (alq_string == "COMP") { + m_alq_type = ALQ_COMP; + } + else if (alq_string == "BEAN") { + m_alq_type = ALQ_BEAN; + } + else if (alq_string == " ") { + m_alq_type = ALQ_UNDEF; + } + else { + m_alq_type = ALQ_INVALID; + throw std::invalid_argument("Invalid ALQ_DEF string"); + } + } + else { + m_alq_type = ALQ_UNDEF; + } + + //Units used for this table + if (itemValid(header, "UNITS")) { + //TODO: Should check that table unit matches rest of deck. + std::string unit_string = header->getItem("UNITS")->getString(0); + if (unit_string == "METRIC") { + } + else if (unit_string == "FIELD") { + } + else if (unit_string == "LAB") { + } + else if (unit_string == "PVT-M") { + } + else { + throw std::invalid_argument("Invalid UNITS string"); + } + } + else { + //Do nothing, table implicitly same unit as rest of deck + } + + //Quantity in the body of the table + if (itemValid(header, "BODY_DEF")) { + std::string body_string = header->getItem("BODY_DEF")->getString(0); + if (body_string == "TEMP") { + throw std::invalid_argument("Invalid BODY_DEF string: TEMP not supported"); + } + else if (body_string == "BHP") { + + } + else { + throw std::invalid_argument("Invalid BODY_DEF string"); + } + } + else { + //Default to BHP + } + + + //Get actual rate / flow values + m_flo_data = (*iter++)->getItem("FLOW_VALUES")->getSIDoubleData(); + + //Get actual tubing head pressure values + m_thp_data = (*iter++)->getItem("THP_VALUES")->getSIDoubleData(); + + //Get actual water fraction values + m_wfr_data = (*iter++)->getItem("WFR_VALUES")->getRawDoubleData(); //FIXME: unit + + //Get actual gas fraction values + m_gfr_data = (*iter++)->getItem("GFR_VALUES")->getRawDoubleData(); //FIXME: unit + + //Get actual gas fraction values + m_alq_data = (*iter++)->getItem("ALQ_VALUES")->getRawDoubleData(); //FIXME: unit + + //Finally, read the actual table itself. + size_t nt = m_thp_data.size(); + size_t nw = m_wfr_data.size(); + size_t ng = m_gfr_data.size(); + size_t na = m_alq_data.size(); + size_t nf = m_flo_data.size(); + extents shape; + shape[0] = nt; + shape[1] = nw; + shape[2] = ng; + shape[3] = na; + shape[4] = nf; + m_data.resize(shape); + + for (; iter!=table->end(); ++iter) { + //Get indices (subtract 1 to get 0-based index) + int t = (*iter)->getItem("THP_INDEX")->getInt(0) - 1; + int w = (*iter)->getItem("WFR_INDEX")->getInt(0) - 1; + int g = (*iter)->getItem("GFR_INDEX")->getInt(0) - 1; + int a = (*iter)->getItem("ALQ_INDEX")->getInt(0) - 1; + + //Rest of values (bottom hole pressure or tubing head temperature) have index of flo value + const std::vector& bhp_tht = (*iter)->getItem("VALUES")->getRawDoubleData(); //FIXME: unit + for (int f=0; f& getFloAxis() { + return m_flo_data; + } + + /** + * Returns the coordinates for the tubing head pressure sample points in the table + * @return Tubing head pressure coordinates + */ + inline const std::vector& getTHPAxis() { + return m_thp_data; + } + + /** + * Returns the coordinates for the water fraction sample points in the table + * @return Water fraction coordinates + */ + inline const std::vector& getWFRAxis() { + return m_wfr_data; + } + + /** + * Returns the coordinates for the gas fraction sample points in the table + * @return Gas fraction coordinates + */ + inline const std::vector& getGFRAxis() { + return m_gfr_data; + } + + /** + * Returns the coordinates for the artificial lift quantity points in the table + * @return Artificial lift quantity coordinates + */ + inline const std::vector& getALQAxis() { + return m_alq_data; + } + + /** + * Returns the data of the table itself. The data is ordered so that + * + * table = getTable(); + * bhp = table[thp_idx][wfr_idx][gfr_idx][alq_idx][flo_idx]; + * + * gives the bottom hole pressure value in the table for the coordinate + * given by + * flo_axis = getFloAxis(); + * thp_axis = getTHPAxis(); + * ... + * + * flo_coord = flo_axis(flo_idx); + * thp_coord = thp_axis(thp_idx); + * ... + */ + inline const array_type& getTable() { + return m_data; + } + +private: + + //"Header" variables + int m_table_num; + double m_datum_depth; + FLO_TYPE m_flo_type; + WFR_TYPE m_wfr_type; + GFR_TYPE m_gfr_type; + ALQ_TYPE m_alq_type; + + //The actual table axes + std::vector m_flo_data; + std::vector m_thp_data; + std::vector m_wfr_data; + std::vector m_gfr_data; + std::vector m_alq_data; + + //The data itself, using the data ordering m_data[thp][wfr][gfr][alq][flo] + array_type m_data; + + + /** + * Helper function that checks if an item exists in a record, and has a + * non-zero size + */ + bool itemValid(DeckRecordConstPtr& record, const char* name) { + if (record->size() == 0) { + return false; + } + else { + DeckItemPtr item; + //TODO: Should we instead here allow the exception to propagate? + try { + item = record->getItem(name); + } + catch (...) { + return false; + } + + if (item->size() > 0) { + return true; + } + else { + return false; + } + } + } + + /** + * Debug function that runs a series of asserts to check for sanity of inputs. + * Called after init to check that everything looks ok. + */ + void check() { + //Table number + assert(m_table_num > 0); + + //Misc types + assert(m_flo_type >= FLO_OIL && m_flo_type < FLO_INVALID); + assert(m_wfr_type >= WFR_WOR && m_wfr_type < WFR_INVALID); + assert(m_gfr_type >= GFR_GOR && m_gfr_type < GFR_INVALID); + assert(m_alq_type >= ALQ_GRAT && m_alq_type < ALQ_INVALID); + + //Data axis size + assert(m_flo_data.size() > 0); + assert(m_thp_data.size() > 0); + assert(m_wfr_data.size() > 0); + assert(m_gfr_data.size() > 0); + assert(m_alq_data.size() > 0); + + //Data axis sorted? + assert(is_sorted(m_flo_data.begin(), m_flo_data.end())); + assert(is_sorted(m_thp_data.begin(), m_thp_data.end())); + assert(is_sorted(m_wfr_data.begin(), m_wfr_data.end())); + assert(is_sorted(m_gfr_data.begin(), m_gfr_data.end())); + assert(is_sorted(m_alq_data.begin(), m_alq_data.end())); + + //Check data size matches axes + assert(m_data.num_dimensions() == 5); + assert(m_data.shape()[0] == m_thp_data.size()); + assert(m_data.shape()[1] == m_wfr_data.size()); + assert(m_data.shape()[2] == m_gfr_data.size()); + assert(m_data.shape()[3] == m_alq_data.size()); + assert(m_data.shape()[4] == m_flo_data.size()); + + //Finally, check that all data is within reasonable ranges, defined to be up-to 1.0e10... + typedef array_type::size_type size_type; + for (size_type t=0; t 1.0e10) { + //TODO: Replace with proper log message + std::cerr << "Too large value encountered in VFPPROD in [" + << t << "," << w << "," << g << "," << a << "," << f << "]=" + << m_data[t][w][g][a][f] << std::endl; + } + } + } + } + } + } + } +}; + +} + + +#endif /* OPM_PARSER_ECLIPSE_ECLIPSESTATE_TABLES_VFPPRODTABLE_HPP_ */ diff --git a/opm/parser/eclipse/EclipseState/Tables/tests/TableTests.cpp b/opm/parser/eclipse/EclipseState/Tables/tests/TableTests.cpp index 290528553..a51f687c7 100644 --- a/opm/parser/eclipse/EclipseState/Tables/tests/TableTests.cpp +++ b/opm/parser/eclipse/EclipseState/Tables/tests/TableTests.cpp @@ -32,6 +32,7 @@ #include #include #include +#include #include @@ -355,3 +356,120 @@ BOOST_AUTO_TEST_CASE(PvtoTable_Tests) { BOOST_CHECK_EQUAL(pvto2OuterTable->getGasSolubilityColumn().front(), 12.0); BOOST_CHECK_EQUAL(pvto2OuterTable->getGasSolubilityColumn().back(), 23.0); } + + +BOOST_AUTO_TEST_CASE(VFPProdTable_Tests) { + const char *deckData = "\ +VFPPROD \n\ +-- Table Depth Rate WFR GFR TAB ALQ UNITS BODY \n\ +-- ----- ----- ----- ----- ----- ----- --- -------- ----- \n\ + 5 32.9 'LIQ' 'WCT' 'GOR' 'THP' ' ' 'METRIC' 'BHP' / \n\ +-- Rate axis \n\ +1 3 5 / \n\ +-- THP axis \n\ +7 11 / \n\ +-- WFR axis \n\ +13 17 / \n\ +-- GFR axis \n\ +19 23 / \n\ +-- ALQ axis \n\ +29 31 / \n\ +-- Table data with THP# WFR# GFR# ALQ# \n\ +1 1 1 1 1.5 2.5 3.5 / \n\ +2 1 1 1 4.5 5.5 6.5 / \n\ +1 2 1 1 7.5 8.5 9.5 / \n\ +2 2 1 1 10.5 11.5 12.5 / \n\ +1 1 2 1 13.5 14.5 15.5 / \n\ +2 1 2 1 16.5 17.5 18.5 / \n\ +1 2 2 1 19.5 20.5 21.5 / \n\ +2 2 2 1 22.5 23.5 24.5 / \n\ +1 1 1 2 25.5 26.5 27.5 / \n\ +2 1 1 2 28.5 29.5 30.5 / \n\ +1 2 1 2 31.5 32.5 33.5 / \n\ +2 2 1 2 34.5 35.5 36.5 / \n\ +1 1 2 2 37.5 38.5 39.5 / \n\ +2 1 2 2 40.5 41.5 42.5 / \n\ +1 2 2 2 43.5 44.5 45.5 / \n\ +2 2 2 2 46.5 47.5 48.5 / \n"; + + Opm::ParserPtr parser(new Opm::Parser); + Opm::DeckConstPtr deck(parser->parseString(deckData)); + Opm::DeckKeywordConstPtr vfpprodKeyword = deck->getKeyword("VFPPROD"); + + BOOST_CHECK_EQUAL(deck->numKeywords("VFPPROD"), 1); + + Opm::VFPProdTable vfpprodTable; + + vfpprodTable.init(vfpprodKeyword); + + BOOST_CHECK_EQUAL(vfpprodTable.getTableNum(), 5); + BOOST_CHECK_EQUAL(vfpprodTable.getDatumDepth(), 32.9); + BOOST_CHECK_EQUAL(vfpprodTable.getFloType(), Opm::VFPProdTable::FLO_LIQ); + BOOST_CHECK_EQUAL(vfpprodTable.getWFRType(), Opm::VFPProdTable::WFR_WCT); + BOOST_CHECK_EQUAL(vfpprodTable.getGFRType(), Opm::VFPProdTable::GFR_GOR); + BOOST_CHECK_EQUAL(vfpprodTable.getALQType(), Opm::VFPProdTable::ALQ_UNDEF); + //TODO: Units + //TODO: Table entries can be either BHP or TEMP => Tubing head temperature. Only BHP supported now + + //Flo axis + { + const std::vector& flo = vfpprodTable.getFloAxis(); + BOOST_REQUIRE_EQUAL(flo.size(), 3); + BOOST_CHECK_EQUAL(flo[0], 1); + BOOST_CHECK_EQUAL(flo[1], 3); + BOOST_CHECK_EQUAL(flo[2], 5); + } + + //THP axis + { + const std::vector& thp = vfpprodTable.getTHPAxis(); + BOOST_REQUIRE_EQUAL(thp.size(), 2); + BOOST_CHECK_EQUAL(thp[0], 700000); + BOOST_CHECK_EQUAL(thp[1], 1100000); + } + + //WFR axis + { + const std::vector& wfr = vfpprodTable.getWFRAxis(); + BOOST_REQUIRE_EQUAL(wfr.size(), 2); + BOOST_CHECK_EQUAL(wfr[0], 13); + BOOST_CHECK_EQUAL(wfr[1], 17); + } + + //GFR axis + { + const std::vector& gfr = vfpprodTable.getGFRAxis(); + BOOST_REQUIRE_EQUAL(gfr.size(), 2); + BOOST_CHECK_EQUAL(gfr[0], 19); + BOOST_CHECK_EQUAL(gfr[1], 23); + } + + //ALQ axis + { + const std::vector& alq = vfpprodTable.getALQAxis(); + BOOST_REQUIRE_EQUAL(alq.size(), 2); + BOOST_CHECK_EQUAL(alq[0], 29); + BOOST_CHECK_EQUAL(alq[1], 31); + } + + //The data itself + { + typedef Opm::VFPProdTable::array_type::size_type size_type; + const Opm::VFPProdTable::array_type& data = vfpprodTable.getTable(); + const size_type* size = data.shape(); + + double index = 0.5; + for (size_type a=0; a