OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge d9aac3dfdd into 9a0ccfd853

Browse Source
This commit is contained in:
Bård Skaflestad 2023-11-24 17:28:53 +01:00 committed by GitHub
commit bf7cb33de1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
12 changed files with 769 additions and 226 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
opm/input/eclipse/Deck/DeckItem.hpp
View File

@ -47,6 +47,7 @@ namespace Opm {
DeckItem( const std::string&, double, const std::vector<Dimension>& active_dim, const std::vector<Dimension>& default_dim);
static DeckItem serializationTestObject();
DeckItem emptyStructuralCopy() const;
const std::string& name() const;
@ -76,6 +77,10 @@ namespace Opm {
template< typename T > const std::vector< T >& getData() const;
const std::vector< double >& getSIDoubleData() const;
const std::vector<value::status>& getValueStatus() const;
const std::vector<Dimension>& getActiveDimensions() const
{
return this->active_dimensions;
}
template< typename T>
void shrink_to_fit();
@ -151,6 +156,7 @@ namespace Opm {
}
void reserve_additionalRawString(std::size_t);
private:
mutable std::vector< double > dval;
std::vector< int > ival;

2
opm/input/eclipse/Deck/DeckKeyword.hpp
View File

@ -51,7 +51,7 @@ namespace Opm {
const std::string& name() const;
void setFixedSize();
const KeywordLocation& location() const;
DeckKeyword emptyStructuralCopy() const;
size_t size() const;
bool empty() const;

2
opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp
View File

@ -39,6 +39,8 @@ namespace Opm {
SimpleTable(TableSchema, const std::string& tableName, const DeckItem& deckItem, const int tableID);
explicit SimpleTable( TableSchema );
virtual ~SimpleTable() = default;
static SimpleTable serializationTestObject();
void addColumns();

103
opm/input/eclipse/EclipseState/Tables/TableContainer.hpp
View File

@ -28,69 +28,77 @@ namespace Opm {
class SimpleTable;
class TableContainer {
/*
The TableContainer class implements a simple map:
{tableNumber , Table}. The main functionality of the
TableContainer class is that the getTable method implements
the Eclipse behavior:
If table N is not implemented - use table N - 1.
The getTable() method will eventually throw an exception if
not even table 0 is there.
Consider the following code:
TableContainer container(10);
std::shared_ptr<TableType> table0 = std::make_shared<TableType>(...);
container.addTable( table0 , 0 )
We create a container with maximum 10 tables, and then we add
one single table at slot 0; then we have:
container.size() == 1
container.hasTable( 0 ) == true
container.hasTable( 9 ) == false
container.hasTable(10 ) == false
container.getTable( 0 ) == container[9] == table0;
container.gteTable(10 ) ==> exception
*/
/// The TableContainer class implements a simple map:
///
/// {tableNumber, Table}
///
/// The main functionality of the TableContainer class is that the
/// getTable() method implements the ECLIPSE behavior:
///
/// If table N is not implemented - use table N - 1.
///
/// The getTable() method will eventually throw an exception if not even
/// table 0 is there.
///
/// Consider the following code:
///
/// TableContainer container(10);
///
/// std::shared_ptr<TableType> table0 = std::make_shared<TableType>(...);
/// container.addTable(table0, 0)
///
/// We create a container with maximum 10 tables, and then we add one
/// single table at slot 0; then we have:
///
/// container.size() == 1
/// container.hasTable( 0 ) == true
/// container.hasTable( 9 ) == false
/// container.hasTable(10 ) == false
///
/// container.getTable( 0 ) == container[9] == table0;
/// container.gteTable(10 ) ==> exception
class TableContainer
{
public:
using TableMap = std::map<size_t, std::shared_ptr<SimpleTable>>;
TableContainer();
explicit TableContainer( size_t maxTables );
explicit TableContainer(size_t maxTables);
static TableContainer serializationTestObject();
bool empty() const;
/*
This is the number of actual tables in the container.
*/
// This is the number of actual tables in the container.
size_t size() const;
size_t max() const;
const TableMap& tables() const;
void addTable(size_t tableNumber , std::shared_ptr<SimpleTable> table);
void addTable(size_t tableNumber, std::shared_ptr<SimpleTable> table);
/*
Observe that the hasTable() method does not invoke the "If
table N is not implemented use table N - 1 behavior.
*/
size_t hasTable(size_t tableNumber) const;
// Observe that the hasTable() method does not invoke the "If table
// N is not implemented use table N - 1 behavior.
bool hasTable(size_t tableNumber) const;
const SimpleTable& getTable(size_t tableNumber) const;
const SimpleTable& operator[](size_t tableNumber) const;
const SimpleTable& operator[](size_t tableNumber) const
{
return this->getTable(tableNumber);
}
template <class TableType>
const TableType& getTable(size_t tableNumber) const {
const SimpleTable &simpleTable = getTable( tableNumber );
const TableType * table = static_cast<const TableType *>( &simpleTable );
return *table;
const TableType& getTable(size_t tableNumber) const
{
// This is, strictly speaking, a downcast so we should prefer
// dynamic_cast<>() instead. However, serializeOp() by
// construction throws away the derived TableType during object
// distribution, keeping only the SimpleTable, so dynamic_cast<>
// will throw a bad_cast exception on ranks other than the I/O
// rank (0). We therefore resort to static_cast<>() here
// instead and hope that the caller specifies the correct
// derived type...
return static_cast<const TableType&>(this->getTable(tableNumber));
}
bool operator==(const TableContainer& data) const;
@ -109,5 +117,4 @@ namespace Opm {
}
#endif
#endif // OPM_TABLE_CONTAINER_HPP

16
src/opm/input/eclipse/Deck/DeckItem.cpp
View File

@ -130,6 +130,22 @@ DeckItem DeckItem::serializationTestObject()
return result;
}
DeckItem DeckItem::emptyStructuralCopy() const
{
auto ret = *this;
ret.dval .clear();
ret.ival .clear();
ret.sval .clear();
ret.rsval.clear();
ret.uval .clear();
ret.value_status.clear();
ret.raw_data = true;
return ret;
}
const std::string& DeckItem::name() const {
return this->item_name;
}

9
src/opm/input/eclipse/Deck/DeckKeyword.cpp
View File

@ -65,6 +65,15 @@ namespace Opm {
return result;
}
DeckKeyword DeckKeyword::emptyStructuralCopy() const
{
auto ret = *this;
ret.m_recordList.clear();
return ret;
}
namespace {
template <typename T>
void add_deckvalue( DeckItem deck_item, DeckRecord& deck_record, const ParserItem& parser_item, const std::vector<DeckValue>& input_record, size_t j) {

221
src/opm/input/eclipse/EclipseState/Tables/SimpleTable.cpp
View File

@ -16,34 +16,41 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <utility>
#include <iostream>
#include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableSchema.hpp>
#include <opm/input/eclipse/Deck/DeckItem.hpp>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
#include <iostream>
#include <fmt/format.h>
namespace Opm {
SimpleTable::SimpleTable( TableSchema schema, const std::string& tableName, const DeckItem& deckItem,
const int tableID) :
m_schema( std::move( schema ) ),
m_jfunc (false)
SimpleTable::SimpleTable(TableSchema schema,
const std::string& tableName,
const DeckItem& deckItem,
const int tableID)
: m_schema(std::move(schema))
, m_jfunc (false)
{
init(tableName, deckItem, tableID );
this->init(tableName, deckItem, tableID);
}
SimpleTable::SimpleTable( TableSchema schema ) :
m_schema( std::move( schema ) ),
m_jfunc (false)
SimpleTable::SimpleTable(TableSchema schema)
: m_schema(std::move(schema))
, m_jfunc (false)
{
addColumns();
this->addColumns();
}
SimpleTable SimpleTable::serializationTestObject()
{
SimpleTable result;
@ -54,44 +61,50 @@ namespace Opm {
return result;
}
void SimpleTable::addRow( const std::vector<double>& row, const std::string& tableName) {
if (row.size() == numColumns()) {
for (size_t colIndex = 0; colIndex < numColumns(); colIndex++) {
auto& col = getColumn( colIndex );
col.addValue( row[colIndex], tableName );
}
} else
void SimpleTable::addRow(const std::vector<double>& row,
const std::string& tableName)
{
const auto ncol = this->numColumns();
if (row.size() != ncol) {
throw std::invalid_argument("Size mismatch");
}
}
void SimpleTable::addColumns() {
for (size_t colIdx = 0; colIdx < m_schema.size(); ++colIdx) {
const auto& schemaColumn = m_schema.getColumn( colIdx );
TableColumn column(schemaColumn); // Some move trickery here ...
m_columns.insert( std::make_pair( schemaColumn.name() , column ));
for (auto colIndex = 0*ncol; colIndex < ncol; ++colIndex) {
this->getColumn(colIndex).addValue(row[colIndex], tableName);
}
}
void SimpleTable::addColumns()
{
const auto ncol = this->m_schema.size();
double SimpleTable::get(const std::string& column , size_t row) const {
const auto& col = getColumn( column );
return col[row];
for (auto colIdx = 0*ncol; colIdx < ncol; ++colIdx) {
const auto& schemaColumn = m_schema.getColumn(colIdx);
this->m_columns.insert(std::make_pair(schemaColumn.name(), TableColumn { schemaColumn }));
}
}
double SimpleTable::get(size_t column , size_t row) const {
const auto& col = getColumn( column );
return col[row];
double SimpleTable::get(const std::string& column, size_t row) const
{
return this->getColumn(column)[row];
}
void SimpleTable::init( const std::string& tableName,
const DeckItem& deckItem,
const int tableID,
double scaling_factor) {
double SimpleTable::get(size_t column, size_t row) const
{
return this->getColumn(column)[row];
}
void SimpleTable::init(const std::string& tableName,
const DeckItem& deckItem,
const int tableID,
const double scaling_factor)
{
this->addColumns();
if ( (deckItem.data_size() % numColumns()) != 0) {
const auto ncol = this->numColumns();
if ((deckItem.data_size() % ncol) != 0) {
throw std::runtime_error {
fmt::format("For table {} with ID {}: "
"Number of input table elements ({}) is "
@ -101,93 +114,117 @@ namespace Opm {
};
}
size_t rows = deckItem.data_size() / numColumns();
for (size_t colIdx = 0; colIdx < numColumns(); ++colIdx) {
auto& column = getColumn( colIdx );
for (size_t rowIdx = 0; rowIdx < rows; rowIdx++) {
size_t deckItemIdx = rowIdx*numColumns() + colIdx;
if (deckItem.defaultApplied(deckItemIdx))
column.addDefault( tableName );
const size_t rows = deckItem.data_size() / ncol;
for (size_t colIdx = 0; colIdx < ncol; ++colIdx) {
auto& column = this->getColumn(colIdx);
for (size_t rowIdx = 0; rowIdx < rows; ++rowIdx) {
const size_t deckItemIdx = rowIdx*ncol + colIdx;
if (deckItem.defaultApplied(deckItemIdx)) {
column.addDefault(tableName);
}
else if (m_jfunc) {
column.addValue( deckItem.getData<double>()[deckItemIdx], tableName );
column.addValue(deckItem.getData<double>()[deckItemIdx], tableName);
}
else if (scaling_factor > 0.0) {
column.addValue(scaling_factor * deckItem.get<double>(deckItemIdx), tableName);
}
else {
if (scaling_factor > 0.0)
column.addValue( scaling_factor * deckItem.get<double>(deckItemIdx), tableName );
else
column.addValue( deckItem.getSIDouble(deckItemIdx), tableName );
column.addValue(deckItem.getSIDouble(deckItemIdx), tableName);
}
}
if (colIdx > 0)
column.applyDefaults(getColumn( 0 ), tableName);
if (colIdx > 0) {
column.applyDefaults(this->getColumn(0), tableName);
}
}
}
size_t SimpleTable::numColumns() const {
size_t SimpleTable::numColumns() const
{
return m_schema.size();
}
size_t SimpleTable::numRows() const {
return getColumn( 0 ).size();
size_t SimpleTable::numRows() const
{
return this->getColumn(0).size();
}
const TableColumn& SimpleTable::getColumn( const std::string& name) const {
if (!this->m_jfunc)
return m_columns.get( name );
const TableColumn& SimpleTable::getColumn(const std::string& name) const
{
if (! this->m_jfunc) {
return this->m_columns.get(name);
}
if (name == "PCOW" || name == "PCOG")
if ((name == "PCOW") || (name == "PCOG")) {
assertJFuncPressure(false); // this will throw since m_jfunc=true
return m_columns.get( name );
}
return this->m_columns.get(name);
}
const TableColumn& SimpleTable::getColumn( size_t columnIndex ) const {
return m_columns.iget( columnIndex );
const TableColumn& SimpleTable::getColumn(size_t columnIndex) const
{
return this->m_columns.iget(columnIndex);
}
TableColumn& SimpleTable::getColumn(const std::string& name)
{
if (! this->m_jfunc) {
return this->m_columns.get(name);
}
TableColumn& SimpleTable::getColumn( const std::string& name) {
if (!this->m_jfunc)
return m_columns.get( name );
if (name == "PCOW" || name == "PCOG")
if ((name == "PCOW") || (name == "PCOG")) {
assertJFuncPressure(false); // this will throw since m_jfunc=true
return m_columns.get( name );
}
return this->m_columns.get(name);
}
TableColumn& SimpleTable::getColumn( size_t columnIndex ) {
return m_columns.iget( columnIndex );
TableColumn& SimpleTable::getColumn(size_t columnIndex)
{
return this->m_columns.iget(columnIndex);
}
bool SimpleTable::hasColumn(const std::string& name) const {
return m_schema.hasColumn( name );
bool SimpleTable::hasColumn(const std::string& name) const
{
return this->m_schema.hasColumn(name);
}
double SimpleTable::evaluate(const std::string& columnName, double xPos) const
{
const auto& argColumn = getColumn( 0 );
const auto& valueColumn = getColumn( columnName );
const auto index = this->getColumn(0).lookup(xPos);
const auto index = argColumn.lookup( xPos );
return valueColumn.eval( index );
return this->getColumn(columnName).eval(index);
}
void SimpleTable::assertJFuncPressure(const bool jf) const {
if (jf == m_jfunc)
void SimpleTable::assertJFuncPressure(const bool jf) const
{
if (jf == this->m_jfunc) {
return;
// if we reach here, wrong values are read from the deck! (JFUNC is used
// incorrectly.) This function writes to std err for now, but will
// after a grace period be rewritten to throw (std::invalid_argument).
if (m_jfunc)
std::cerr << "Developer warning: Pressure column is read with JFUNC in deck." << std::endl;
else
std::cerr << "Developer warning: Raw values from JFUNC column is read, but JFUNC not provided in deck." << std::endl;
}
// if we reach here, wrong values are read from the deck! (JFUNC is
// used incorrectly.) This function writes to std err for now, but
// will after a grace period be rewritten to throw
// (std::invalid_argument).
if (m_jfunc) {
std::cerr << "Developer warning: Pressure column "
"is read with JFUNC in deck.\n";
}
else {
std::cerr << "Developer warning: Raw values from JFUNC column "
"is read, but JFUNC not provided in deck.\n";
}
}
bool SimpleTable::operator==(const SimpleTable& data) const {
return this->m_schema == data.m_schema &&
this->m_columns == data.m_columns &&
this->m_jfunc == data.m_jfunc;
bool SimpleTable::operator==(const SimpleTable& data) const
{
return (this->m_schema == data.m_schema)
&& (this->m_columns == data.m_columns)
&& (this->m_jfunc == data.m_jfunc)
;
}
}

111
src/opm/input/eclipse/EclipseState/Tables/TableContainer.cpp
View File

@ -17,22 +17,24 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableContainer.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
#include <string>
#include <utility>
#include <stddef.h>
namespace Opm {
TableContainer::TableContainer() :
m_maxTables(0)
{
}
TableContainer::TableContainer()
: m_maxTables(0)
{}
TableContainer::TableContainer(size_t maxTables) :
m_maxTables(maxTables)
{
}
TableContainer::TableContainer(size_t maxTables)
: m_maxTables(maxTables)
{}
TableContainer TableContainer::serializationTestObject()
{
@ -44,77 +46,86 @@ namespace Opm {
return result;
}
bool TableContainer::empty() const {
bool TableContainer::empty() const
{
return m_tables.empty();
}
size_t TableContainer::size() const {
size_t TableContainer::size() const
{
return m_tables.size();
}
size_t TableContainer::max() const {
size_t TableContainer::max() const
{
return m_maxTables;
}
const TableContainer::TableMap& TableContainer::tables() const {
const TableContainer::TableMap& TableContainer::tables() const
{
return m_tables;
}
size_t TableContainer::hasTable(size_t tableNumber) const {
if (m_tables.find( tableNumber ) == m_tables.end())
return false;
else
return true;
bool TableContainer::hasTable(size_t tableNumber) const
{
return this->m_tables.find(tableNumber)
!= this->m_tables.end();
}
const SimpleTable& TableContainer::getTable(size_t tableNumber) const {
if (tableNumber >= m_maxTables)
const SimpleTable& TableContainer::getTable(size_t tableNumber) const
{
if (tableNumber >= m_maxTables) {
throw std::invalid_argument("TableContainer - invalid tableNumber");
}
if (hasTable(tableNumber)) {
auto pair = m_tables.find( tableNumber );
return *(pair->second.get());
} else {
if (tableNumber > 0)
return getTable(tableNumber -1);
else
throw std::invalid_argument("TableContainer does not have any table in the range 0..." + std::to_string( tableNumber ));
if (auto pairPos = m_tables.find(tableNumber); pairPos != m_tables.end()) {
return *pairPos->second;
}
else if (tableNumber > 0) {
return getTable(tableNumber - 1);
}
else {
throw std::invalid_argument {
"TableContainer does not have any table in "
"the range 0..." + std::to_string(tableNumber)
};
}
}
void TableContainer::addTable(size_t tableNumber, std::shared_ptr<SimpleTable> table)
{
if (tableNumber >= m_maxTables) {
throw std::invalid_argument {
"TableContainer has max: " + std::to_string(m_maxTables) +
" tables. Table number: " + std::to_string(tableNumber) + " illegal."
};
}
const SimpleTable& TableContainer::operator[](size_t tableNumber) const {
return getTable(tableNumber);
m_tables[tableNumber] = std::move(table);
}
void TableContainer::addTable(size_t tableNumber , std::shared_ptr<SimpleTable> table) {
if (tableNumber >= m_maxTables)
throw std::invalid_argument("TableContainer has max: " + std::to_string( m_maxTables ) + " tables. Table number: " + std::to_string( tableNumber ) + " illegal.");
m_tables[tableNumber] = table;
}
bool TableContainer::operator==(const TableContainer& data) const {
if (this->max() != data.max())
bool TableContainer::operator==(const TableContainer& data) const
{
if (this->max() != data.max()) {
return false;
if (this->size() != data.size())
}
if (this->size() != data.size()) {
return false;
}
for (const auto& it : m_tables) {
auto it2 = data.m_tables.find(it.first);
if (it2 == data.m_tables.end())
if (it2 == data.m_tables.end()) {
return false;
if (!(*it.second == *it2->second))
}
if (! (*it.second == *it2->second)) {
return false;
}
}
return true;
}
}
} // namespace Opm

21
src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp
View File

@ -355,8 +355,9 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
}
void TableManager::addTables( const std::string& tableName , size_t numTables) {
m_simpleTables.emplace(std::make_pair(tableName , TableContainer( numTables )));
void TableManager::addTables(const std::string& tableName, size_t numTables)
{
this->m_simpleTables.try_emplace(tableName, numTables);
}
@ -1555,24 +1556,26 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
void TableManager::initSimpleTableContainer(const Deck& deck,
const std::string& keywordName,
const std::string& tableName,
size_t numTables) {
if (!deck.hasKeyword(keywordName))
size_t numTables)
{
if (!deck.hasKeyword(keywordName)) {
return; // the table is not featured by the deck...
auto& container = forceGetTables(tableName , numTables);
}
if (deck.count(keywordName) > 1) {
complainAboutAmbiguousKeyword(deck, keywordName);
return;
}
auto& container = forceGetTables(tableName, numTables);
auto lastComplete = 0 * numTables;
const auto& tableKeyword = deck[keywordName].back();
for (size_t tableIdx = 0; tableIdx < tableKeyword.size(); ++tableIdx) {
const auto& dataItem = tableKeyword.getRecord( tableIdx ).getItem("DATA");
const auto& dataItem = tableKeyword.getRecord(tableIdx).getItem("DATA");
if (dataItem.data_size() > 0) {
std::shared_ptr<TableType> table = std::make_shared<TableType>( dataItem, tableIdx );
container.addTable( tableIdx , table );
container.addTable(tableIdx, std::make_shared<TableType>(dataItem, tableIdx));
lastComplete = tableIdx;
}
else if (tableIdx > static_cast<size_t>(0)) {

489
src/opm/input/eclipse/EclipseState/Tables/Tables.cpp
View File

@ -15,23 +15,13 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
*/
#include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
#include <opm/input/eclipse/Deck/Deck.hpp>
#include <opm/input/eclipse/Deck/DeckItem.hpp>
#include <opm/input/eclipse/Deck/DeckKeyword.hpp>
#include <opm/input/eclipse/Deck/DeckRecord.hpp>
#include <opm/input/eclipse/EclipseState/Tables/ColumnSchema.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableColumn.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableSchema.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/C.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/D.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/G.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/P.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/R.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/V.hpp>
#include <opm/input/eclipse/Units/Dimension.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/input/eclipse/EclipseState/Tables/AqutabTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/EnkrvdTable.hpp>
@ -40,6 +30,7 @@
#include <opm/input/eclipse/EclipseState/Tables/FoamadsTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/FoammobTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/GasvisctTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/GsfTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/ImkrvdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/ImptvdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/MiscTable.hpp>
@ -91,23 +82,397 @@
#include <opm/input/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TracerVdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/WatvisctTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/GsfTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/WsfTable.hpp>
#include <opm/input/eclipse/Units/Dimension.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/common/utility/OpmInputError.hpp>
#include <opm/input/eclipse/Deck/Deck.hpp>
#include <opm/input/eclipse/Deck/DeckItem.hpp>
#include <opm/input/eclipse/Deck/DeckKeyword.hpp>
#include <opm/input/eclipse/Deck/DeckRecord.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/C.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/D.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/G.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/P.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/R.hpp>
#include <opm/input/eclipse/Parser/ParserKeywords/V.hpp>
#include <algorithm>
#include <array>
#include <cstddef>
#include <functional>
#include <initializer_list>
#include <stdexcept>
#include <string_view>
#include <utility>
#include <vector>
#include <stddef.h>
#include <fmt/format.h>
namespace {
/// Simple query interface for extracting tabulated values of saturated
/// gas.
///
/// In particular, supports querying the tabulated gas pressure and its
/// associated formation volume factor, viscosity and vaporised oil
/// concentration ("Rv") values.
struct GasPropertyTableInterface
{
/// Retrieve gas pressure at saturated conditions
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
virtual double pressure(const std::size_t row) const = 0;
/// Retrieve formation volume factor for gas at saturated
/// conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
virtual double fvf(const std::size_t row) const = 0;
/// Retrieve phase viscosity for gas at saturated conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
virtual double viscosity(const std::size_t row) const = 0;
/// Retrieve vaporised oil concentration ("Rv") for gas at saturated
/// conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
virtual double vaporisedOil(const std::size_t row) const = 0;
};
/// Padding for gas property tables at low pressure
class LowPressureTablePadding
{
public:
/// Constructor
///
/// \param[in] Tabulated gas properties at saturated conditions
explicit LowPressureTablePadding(const GasPropertyTableInterface& prop);
/// Whether or not input table needs padding at low pressures
bool inputNeedsPadding() const { return this->needPadding_; }
/// Low pressure padding rows for input table. Not needed unless
/// inputNeedsPadding() returns \c true.
Opm::SimpleTable padding() const;
private:
/// Gas pressure values in rows zero and one of input table.
std::array<double, 2> p_{};
/// Interpolated gas property values at "limiting" pressure
struct {
/// Limiting pressure
double p {0.0};
/// Vaporised oil concentration at limiting pressure
double rv {0.0};
/// Formation volume factor at limiting pressure
double fvf {0.0};
/// Gas viscosity at limiting pressure
double mu {0.0};
} limit_{};
/// Whether or not input table needs padding at low pressure values.
bool needPadding_{false};
};
LowPressureTablePadding::LowPressureTablePadding(const GasPropertyTableInterface& prop)
: p_{ prop.pressure(0), prop.pressure(1) }
{
auto linInterp = [](const std::array<double, 2>& xi,
const std::array<double, 2>& yi,
const double x)
{
const auto t = (x - xi[0]) / (xi[1] - xi[0]);
return (1.0 - t)*yi[0] + t*yi[1];
};
const auto rv = std::array {
prop.vaporisedOil(0),
prop.vaporisedOil(1),
};
const auto b = std::array {
1.0 / prop.fvf(0), // 1 / B(0)
1.0 / prop.fvf(1), // 1 / B(1)
};
const auto recipBmu = std::array {
b[0] / prop.viscosity(0), // 1 / (B(0) * mu(0))
b[1] / prop.viscosity(1), // 1 / (B(1) * mu(1))
};
this->limit_.fvf = std::max(10.0 / b[0], 1.0);
this->limit_.p = linInterp(b, this->p_, 1.0 / this->limit_.fvf);
this->limit_.rv = linInterp(this->p_, rv, this->limit_.p);
this->limit_.mu = 1.0 /
(this->limit_.fvf * linInterp(this->p_, recipBmu, this->limit_.p));
const auto p0 = 1.0*Opm::unit::barsa;
this->needPadding_ = (this->p_[0] > p0) &&
(linInterp(this->p_, b, p0) < 1.0 / this->limit_.fvf);
}
Opm::SimpleTable LowPressureTablePadding::padding() const
{
auto padSchema = Opm::TableSchema{};
padSchema.addColumn(Opm::ColumnSchema("PG" , Opm::Table::STRICTLY_INCREASING, Opm::Table::DEFAULT_NONE));
padSchema.addColumn(Opm::ColumnSchema("RV" , Opm::Table::RANDOM, Opm::Table::DEFAULT_NONE));
padSchema.addColumn(Opm::ColumnSchema("BG" , Opm::Table::RANDOM, Opm::Table::DEFAULT_LINEAR));
padSchema.addColumn(Opm::ColumnSchema("MUG", Opm::Table::RANDOM, Opm::Table::DEFAULT_LINEAR));
auto padTable = Opm::SimpleTable { std::move(padSchema) };
const auto p0 = 1.0*Opm::unit::barsa;
if (this->limit_.p < this->p_[0]) {
if (p0 < this->limit_.p) {
padTable.addRow({ p0, this->limit_.rv, 1.1*this->limit_.fvf, this->limit_.mu }, "PAD");
}
padTable.addRow({ this->limit_.p, this->limit_.rv, this->limit_.fvf, this->limit_.mu }, "PAD");
}
return padTable;
}
// -------------------------------------------------------------------------
/// Gas property query interface implementation for dry gas
class DryGasTable : public GasPropertyTableInterface
{
public:
/// Constructor
///
/// \param[in] pvdg Tabulated gas property values for dry gas
explicit DryGasTable(const Opm::PvdgTable& pvdg)
: pvdg_ { pvdg }
{}
/// Retrieve gas pressure at saturated conditions
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
double pressure(const std::size_t row) const override;
/// Retrieve formation volume factor for gas at saturated
/// conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
double fvf(const std::size_t row) const override;
/// Retrieve phase viscosity for gas at saturated conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
double viscosity(const std::size_t row) const override;
/// Retrieve vaporised oil concentration ("Rv") for gas at saturated
/// conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case. Ignored
/// for dry gas, as the vaporised oil concentration is always zero
/// in this case.
double vaporisedOil([[maybe_unused]] const std::size_t row) const override
{
return 0.0;
}
private:
/// Underlying property table.
std::reference_wrapper<const Opm::PvdgTable> pvdg_;
};
double DryGasTable::pressure(const std::size_t row) const
{
return this->pvdg_.get().getPressureColumn()[row];
}
double DryGasTable::fvf(const std::size_t row) const
{
return this->pvdg_.get().getFormationFactorColumn()[row];
}
double DryGasTable::viscosity(const std::size_t row) const
{
return this->pvdg_.get().getViscosityColumn()[row];
}
// -------------------------------------------------------------------------
/// Gas property query interface implementation for wet gas
class WetGasTable : public GasPropertyTableInterface
{
public:
/// Constructor
///
/// \param[in] pvtg Tabulated gas property values for wet gas
explicit WetGasTable(const Opm::PvtgTable& pvtg)
: satTable_ { pvtg.getSaturatedTable() }
{}
/// Retrieve gas pressure at saturated conditions
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
double pressure(const std::size_t row) const override;
/// Retrieve formation volume factor for gas at saturated
/// conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
double fvf(const std::size_t row) const override;
/// Retrieve phase viscosity for gas at saturated conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
double viscosity(const std::size_t row) const override;
/// Retrieve vaporised oil concentration ("Rv") for gas at saturated
/// conditions.
///
/// \param[in] row Row index in saturated gas property table.
/// Typically 0 or 1 in the intended/primary use case.
double vaporisedOil(const std::size_t row) const override;
private:
/// Underlying property table for wet gas at saturated conditions.
std::reference_wrapper<const Opm::SimpleTable> satTable_;
};
double WetGasTable::pressure(const std::size_t row) const
{
return this->satTable_.get().getColumn(0)[row];
}
double WetGasTable::fvf(const std::size_t row) const
{
return this->satTable_.get().getColumn(2)[row];
}
double WetGasTable::viscosity(const std::size_t row) const
{
return this->satTable_.get().getColumn(3)[row];
}
double WetGasTable::vaporisedOil(const std::size_t row) const
{
return this->satTable_.get().getColumn(1)[row];
}
} // Anonymous namespace
namespace Opm
{
namespace {
DeckItem createPvtgItemZero(const DeckKeyword& pvtgTableInput)
{
return pvtgTableInput[0].getItem(0).emptyStructuralCopy();
}
DeckItem createPvtgItemOne(const DeckKeyword& pvtgTableInput)
{
return pvtgTableInput[0].getItem(1).emptyStructuralCopy();
}
DeckKeyword paddedPVTGTable(const SimpleTable& padding,
const DeckKeyword& pvtgTableInput,
const PvtgTable& pvtgTable)
{
auto paddedTable = pvtgTableInput.emptyStructuralCopy();
// Padding
{
const auto nrow = padding.numRows();
const auto ncol = padding.numColumns();
for (auto row = 0*nrow; row < nrow; ++row) {
auto recordItems = std::vector {
createPvtgItemZero(pvtgTableInput),
createPvtgItemOne (pvtgTableInput)
};
{
auto& pg = recordItems.front();
const auto& dim = pg.getActiveDimensions();
pg.push_back(dim.front().convertSiToRaw(padding.get(0, row)));
}
{
auto& rest = recordItems.back();
const auto& dim = rest.getActiveDimensions();
auto n = 0*dim.size();
for (auto col = 1 + 0*ncol; col < ncol; ++col, n = (n + 1) % dim.size()) {
rest.push_back(dim[n].convertSiToRaw(padding.get(col, row)));
}
}
paddedTable.addRecord(DeckRecord { std::move(recordItems) });
}
}
// Original input table
for (auto pgIx = 0*pvtgTable.size(); pgIx < pvtgTable.size(); ++pgIx) {
auto recordItems = std::vector {
createPvtgItemZero(pvtgTableInput),
createPvtgItemOne (pvtgTableInput)
};
{
auto& pg = recordItems.front();
const auto& dim = pg.getActiveDimensions();
pg.push_back(dim.front().convertSiToRaw(pvtgTable.getArgValue(pgIx)));
}
{
auto& rest = recordItems.back();
const auto& dim = rest.getActiveDimensions();
const auto& pgTable = pvtgTable.getUnderSaturatedTable(pgIx);
const auto nrow = pgTable.numRows();
const auto ncol = pgTable.numColumns();
for (auto row = 0*nrow; row < nrow; ++row) {
for (auto col = 0*ncol; col < ncol; ++col) {
rest.push_back(dim[col].convertSiToRaw(pgTable.get(col, row)));
}
}
}
paddedTable.addRecord(DeckRecord { std::move(recordItems) });
}
// Resulting padded table holds values for just a single table/PVT
// region, even if pvtgTableInput holds tables for multiple PVT
// regions.
return paddedTable;
}
} // Anonymous namespace
PvtgTable::PvtgTable(const DeckKeyword& keyword, size_t tableIdx)
: PvtxTable("P")
{
@ -120,7 +485,34 @@ PvtgTable::PvtgTable(const DeckKeyword& keyword, size_t tableIdx)
m_saturatedSchema.addColumn(ColumnSchema("BG", Table::RANDOM, Table::DEFAULT_LINEAR));
m_saturatedSchema.addColumn(ColumnSchema("MUG", Table::RANDOM, Table::DEFAULT_LINEAR));
// Run full table initialisation first. The downside of this is that we
// will end up throwing away and redoing that work if the table needs
// padding at low pressure values. On the other hand, the full table
// initialisation procedure also checks consistency and replaces any
// defaulted values which means we don't have to have special logic in
// place to handle those complexities if the table does need padding.
PvtxTable::init(keyword, tableIdx);
if (const auto tablePadding = LowPressureTablePadding { WetGasTable { *this } };
tablePadding.inputNeedsPadding())
{
// Note: The padded PVTG table holds values for a single table/PVT
// region, even if 'keyword' holds tables for multiple PVT regions.
// We therefore unconditionally pass '0' as the 'tableIdx' in this
// case. Moreover, PvtxTable::init() expects that both the outer
// column and the array of undersaturated tables are empty, so clear
// those here, once we've used their contents to form the padding
// table.
const auto paddedTable =
paddedPVTGTable(tablePadding.padding(), keyword, *this);
this->m_outerColumn = TableColumn { this->m_outerColumnSchema };
this->m_underSaturatedTables.clear();
PvtxTable::init(paddedTable, 0);
}
}
PvtgTable
@ -601,15 +993,78 @@ Sof3Table::getKrogColumn() const
return SimpleTable::getColumn(2);
}
namespace {
DeckItem paddedPVDGTable(const SimpleTable& padding,
const DeckItem& pvdgTableInput,
const PvdgTable& pvdgTable)
{
auto paddedTable = pvdgTableInput.emptyStructuralCopy();
auto addTableValue = [&paddedTable,
&dim = pvdgTableInput.getActiveDimensions(),
n = std::size_t{0}]
(const double value) mutable
{
paddedTable.push_back(dim[n].convertSiToRaw(value));
n = (n + 1) % dim.size();
};
// Padding
{
const auto columns = std::array {
&padding.getColumn("PG"),
&padding.getColumn("BG"),
&padding.getColumn("MUG"),
};
const auto nrow = padding.numRows();
for (auto row = 0*nrow; row < nrow; ++row) {
for (auto* column : columns) {
addTableValue((*column)[row]);
}
}
}
// Original table
{
const auto nrow = pvdgTable.numRows();
const auto ncol = pvdgTable.numColumns();
for (auto row = 0*nrow; row < nrow; ++row) {
for (auto col = 0*ncol; col < ncol; ++col) {
// Yes, argument order is (col, row)
addTableValue(pvdgTable.get(col, row));
}
}
}
return paddedTable;
}
} // Anonymous namespace
PvdgTable::PvdgTable(const DeckItem& item, const int tableID)
{
m_schema.addColumn(ColumnSchema("P", Table::STRICTLY_INCREASING, Table::DEFAULT_NONE));
m_schema.addColumn(ColumnSchema("BG", Table::STRICTLY_DECREASING, Table::DEFAULT_LINEAR));
m_schema.addColumn(ColumnSchema("MUG", Table::INCREASING, Table::DEFAULT_LINEAR));
SimpleTable::init("PVDG", item, tableID);
}
// Run full table initialisation first. The downside of this is that we
// will end up throwing away and redoing that work if the table needs
// padding at low pressure values. On the other hand, the full table
// initialisation procedure also checks consistency and replaces any
// defaulted values which means we don't have to have special logic in
// place to handle those complexities if the table does need padding.
const auto tableName = std::string { "PVDG" };
SimpleTable::init(tableName, item, tableID);
if (const auto tablePadding = LowPressureTablePadding { DryGasTable { *this } };
tablePadding.inputNeedsPadding())
{
SimpleTable::init(tableName, paddedPVDGTable(tablePadding.padding(), item, *this), tableID);
}
}
const TableColumn&
PvdgTable::getPressureColumn() const

9
tests/material/test_eclblackoilpvt.cpp
View File

@ -233,14 +233,11 @@ void ensurePvtApi(const OilPvt& oilPvt, const GasPvt& gasPvt, const WaterPvt& wa
struct Fixture {
Fixture()
: python(std::make_shared<Opm::Python>())
, deck(Opm::Parser().parseString(deckString1))
: deck(Opm::Parser().parseString(deckString1))
, eclState(deck)
, schedule(deck, eclState, python)
{
}
, schedule(deck, eclState, std::make_shared<Opm::Python>())
{}
std::shared_ptr<Opm::Python> python;
Opm::Deck deck;
Opm::EclipseState eclState;
Opm::Schedule schedule;

6
tests/msim/uda.include
View File

@ -1,4 +1,4 @@
std::string uda_deck = R"(
const std::string uda_deck { R"(
-- This reservoir simulation deck is made available under the Open Database
-- License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in
-- individual contents of the database are licensed under the Database Contents
@ -35,6 +35,7 @@ OIL
GAS
WATER
DISGAS
FIELD
START
1 'DEC' 2014 /
@ -380,6 +381,5 @@ DATES
1 'DEC' 2016 /
/
END
)";
)" };