correct the GasvisctTable class

it turns out that the number of columns of this class depends on the
number of specified components. this lead to all kinds of fun like
manual unit conversion, etc.
This commit is contained in:
Andreas Lauser
2015-01-30 14:24:01 +01:00
parent 0a477ed18a
commit dcc36df62e
6 changed files with 123 additions and 17 deletions

View File

@@ -261,7 +261,6 @@ namespace Opm {
void EclipseState::initTables(DeckConstPtr deck, LoggerPtr logger) {
initSimpleTables(deck, logger, "ENKRVD", m_enkrvdTables);
initSimpleTables(deck, logger, "ENPTVD", m_enptvdTables);
initSimpleTables(deck, logger, "GASVISCT", m_gasvisctTables);
initSimpleTables(deck, logger, "IMKRVD", m_imkrvdTables);
initSimpleTables(deck, logger, "IMPTVD", m_imptvdTables);
initSimpleTables(deck, logger, "OILVISCT", m_oilvisctTables);
@@ -279,6 +278,10 @@ namespace Opm {
initSimpleTables(deck, logger, "SWFN", m_swfnTables);
initSimpleTables(deck, logger, "WATVISCT", m_watvisctTables);
// the number of columns of the GASVSISCT tables depends on the value of the
// COMPS keyword...
initGasvisctTables(deck, logger, "GASVISCT", m_gasvisctTables);
// the ROCKTAB table comes with additional fun because the number of columns
//depends on the presence of the RKTRMDIR keyword...
initRocktabTables(deck, logger);
@@ -470,6 +473,39 @@ namespace Opm {
}
}
void EclipseState::initGasvisctTables(DeckConstPtr deck,
LoggerPtr logger,
const std::string& keywordName,
std::vector<GasvisctTable>& tableVector) {
if (!deck->hasKeyword(keywordName))
return; // the table is not featured by the deck...
if (deck->numKeywords(keywordName) > 1) {
complainAboutAmbiguousKeyword(deck, logger, keywordName);
return;
}
const auto& tableKeyword = deck->getKeyword(keywordName);
for (size_t tableIdx = 0; tableIdx < tableKeyword->size(); ++tableIdx) {
if (tableKeyword->getRecord(tableIdx)->getItem(0)->size() == 0) {
// for simple tables, an empty record indicates that the previous table
// should be copied...
if (tableIdx == 0) {
logger->addError(tableKeyword->getFileName(),
tableKeyword->getLineNumber(),
"The first table for keyword "+keywordName+
" must be explicitly defined! Ignoring keyword");
return;
}
tableVector.push_back(tableVector.back());
continue;
}
tableVector.push_back(GasvisctTable());
tableVector[tableIdx].init(deck, tableKeyword, tableIdx);
}
}
bool EclipseState::supportsGridProperty(const std::string& keyword, int enabledTypes) const {
bool result = false;
if (enabledTypes & IntProperties)

View File

@@ -191,6 +191,10 @@ namespace Opm {
}
void initRocktabTables(DeckConstPtr deck, LoggerPtr logger);
void initGasvisctTables(DeckConstPtr deck,
LoggerPtr logger,
const std::string& keywordName,
std::vector<GasvisctTable>& tableVector);
void setMULTFLT(std::shared_ptr<const Section> section, LoggerPtr logger) const;
void initMULTREGT(DeckConstPtr deck, LoggerPtr logger);

View File

@@ -35,22 +35,72 @@ namespace Opm {
* \brief Read the GASVISCT keyword and provide some convenience
* methods for it.
*/
void init(Opm::DeckKeywordConstPtr keyword, int recordIdx)
void init(Opm::DeckConstPtr deck, Opm::DeckKeywordConstPtr keyword, int recordIdx)
{
ParentType::init(keyword,
std::vector<std::string>{
"Temperature",
"Viscosity",
"Unknown" // the RM mentions a third in the example but nowhere else. WTF?
},
recordIdx,
/*firstEntityOffset=*/0);
int numComponents = deck->getKeyword("COMPS")->getRecord(0)->getItem(0)->getInt(0);
auto temperatureDimension = deck->getActiveUnitSystem()->getDimension("Temperature");
auto viscosityDimension = deck->getActiveUnitSystem()->getDimension("Viscosity");
// create the columns: temperature plus one viscosity column per component
std::vector<std::string> columnNames;
columnNames.push_back("Temperature");
for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
columnNames.push_back("Viscosity" + std::to_string(compIdx));
ParentType::createColumns(columnNames);
// extract the actual data from the deck
Opm::DeckRecordConstPtr deckRecord =
keyword->getRecord(recordIdx);
size_t numFlatItems = getNumFlatItems(deckRecord);
if ( numFlatItems % numColumns() != 0)
throw std::runtime_error("Number of columns in the data file is inconsistent "
"with the expected number for keyword GASVISCT");
for (size_t rowIdx = 0; rowIdx*numColumns() < numFlatItems; ++rowIdx) {
// add the current temperature
int deckItemIdx = rowIdx*numColumns();
bool isDefaulted = this->getFlatIsDefaulted(deckRecord, deckItemIdx);
this->m_valueDefaulted[0].push_back(isDefaulted);
if (!isDefaulted) {
double T = this->getFlatRawDoubleData(deckRecord, deckItemIdx);
this->m_columns[0].push_back(temperatureDimension->convertRawToSi(T));
}
// deal with the component viscosities
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
deckItemIdx = rowIdx*numColumns() + compIdx + 1;
size_t columnIdx = compIdx + 1;
isDefaulted = this->getFlatIsDefaulted(deckRecord, deckItemIdx);
this->m_valueDefaulted[columnIdx].push_back(isDefaulted);
if (!isDefaulted) {
double mu = this->getFlatRawDoubleData(deckRecord, deckItemIdx);
this->m_columns[columnIdx].push_back(viscosityDimension->convertRawToSi(mu));
}
}
}
// make sure that the columns agree with the keyword specification of the
// reference manual. (actually, the documentation does not say anyting about
// whether items of these columns are defaultable or not, so we assume here
// that they are not.)
ParentType::checkNonDefaultable("Temperature");
ParentType::checkMonotonic("Temperature", /*isAscending=*/true);
ParentType::checkNonDefaultable("Viscosity");
ParentType::checkMonotonic("Viscosity", /*isAscending=*/true, /*strictlyMonotonic=*/false);
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
std::string columnName = "Viscosity" + std::to_string(compIdx);
ParentType::checkNonDefaultable(columnName);
ParentType::checkMonotonic(columnName,
/*isAscending=*/true,
/*strictlyMonotonic=*/false);
}
}
public:
@@ -62,8 +112,8 @@ namespace Opm {
const std::vector<double> &getTemperatureColumn() const
{ return ParentType::getColumn(0); }
const std::vector<double> &getGasViscosityColumn() const
{ return ParentType::getColumn(1); }
const std::vector<double> &getGasViscosityColumn(size_t compIdx) const
{ return ParentType::getColumn(1 + compIdx); }
};
}

View File

@@ -243,6 +243,20 @@ size_t SingleRecordTable::getNumFlatItems(Opm::DeckRecordConstPtr deckRecord) co
return result;
}
double SingleRecordTable::getFlatRawDoubleData(Opm::DeckRecordConstPtr deckRecord, size_t flatItemIdx) const
{
size_t itemFirstFlatIdx = 0;
for (unsigned i = 0; i < deckRecord->size(); ++ i) {
Opm::DeckItemConstPtr item = deckRecord->getItem(i);
if (itemFirstFlatIdx + item->size() > flatItemIdx)
return item->getRawDouble(flatItemIdx - itemFirstFlatIdx);
else
itemFirstFlatIdx += item->size();
}
throw std::range_error("Tried to access out-of-range flat item");
}
double SingleRecordTable::getFlatSiDoubleData(Opm::DeckRecordConstPtr deckRecord, size_t flatItemIdx) const
{
size_t itemFirstFlatIdx = 0;

View File

@@ -75,6 +75,7 @@ namespace Opm {
* X coordinate.
*/
double evaluate(const std::string& columnName, double xPos) const;
protected:
void checkNonDefaultable(const std::string& columnName);
void checkMonotonic(const std::string& columnName,
@@ -84,6 +85,7 @@ namespace Opm {
void applyDefaultsLinear(const std::string& columnName);
void createColumns(const std::vector<std::string> &columnNames);
size_t getNumFlatItems(Opm::DeckRecordConstPtr deckRecord) const;
double getFlatRawDoubleData(Opm::DeckRecordConstPtr deckRecord, size_t flatItemIdx) const;
double getFlatSiDoubleData(Opm::DeckRecordConstPtr deckRecord, size_t flatItemIdx) const;
bool getFlatIsDefaulted(Opm::DeckRecordConstPtr deckRecord, size_t flatItemIdx) const;