riTRANSXYZ: nnc calculations in place

Making this feaure complete.
Had to add the concept of several results sets to the nnc data
This commit is contained in:
Jacob Støren 2014-08-26 12:07:08 +02:00
parent 56b1f78f2f
commit 670ca4ced5
10 changed files with 255 additions and 33 deletions

View File

@ -479,13 +479,14 @@ void RifReaderEclipseOutput::transferNNCData( const ecl_grid_type * mainEclGrid
//cvf::Trace::show("Reading NNC. Count: " + cvf::String(numNNC));
mainGrid->nncData()->connections().resize(numNNC);
std::vector<double>& transmissibilityValues = mainGrid->nncData()->makeConnectionScalarResult(cvf::UNDEFINED_SIZE_T);
for (int nIdx = 0; nIdx < numNNC; ++nIdx)
{
RigGridBase* grid1 = mainGrid->gridByIndex(eclNNCData[nIdx].grid_nr1);
mainGrid->nncData()->connections()[nIdx].m_c1GlobIdx = grid1->reservoirCellIndex(eclNNCData[nIdx].global_index1);
RigGridBase* grid2 = mainGrid->gridByIndex(eclNNCData[nIdx].grid_nr2);
mainGrid->nncData()->connections()[nIdx].m_c2GlobIdx = grid2->reservoirCellIndex(eclNNCData[nIdx].global_index2);
mainGrid->nncData()->connections()[nIdx].m_transmissibility = eclNNCData[nIdx].trans;
transmissibilityValues[nIdx] = eclNNCData[nIdx].trans;
}

View File

@ -647,11 +647,15 @@ void RivFaultPartMgr::updateNNCColors(RimResultSlot* cellResultSlot)
{
if (m_NNCFaces.isNull()) return;
if (cellResultSlot && cellResultSlot->resultVariable() == RimDefines::combinedTransmissibilityResultName())
if (cellResultSlot
&& ( cellResultSlot->resultVariable() == RimDefines::combinedTransmissibilityResultName()
|| cellResultSlot->resultVariable() == RimDefines::combinedRiTransResultName()))
{
size_t scalarSetIndex = cellResultSlot->gridScalarIndex();
const cvf::ScalarMapper* mapper = cellResultSlot->legendConfig()->scalarMapper();
m_NNCGenerator->textureCoordinates(m_NNCTextureCoords.p(), mapper);
m_NNCGenerator->textureCoordinates(m_NNCTextureCoords.p(), mapper, scalarSetIndex);
cvf::ref<cvf::Effect> nncEffect;

View File

@ -144,12 +144,19 @@ void RivNNCGeometryGenerator::computeArrays()
/// Calculates the texture coordinates in a "nearly" one dimensional texture.
/// Undefined values are coded with a y-texture coordinate value of 1.0 instead of the normal 0.5
//--------------------------------------------------------------------------------------------------
void RivNNCGeometryGenerator::textureCoordinates(cvf::Vec2fArray* textureCoords, const cvf::ScalarMapper* mapper) const
void RivNNCGeometryGenerator::textureCoordinates(cvf::Vec2fArray* textureCoords, const cvf::ScalarMapper* mapper, size_t scalarResultIndex) const
{
size_t numVertices = m_vertices->size();
textureCoords->resize(numVertices);
cvf::Vec2f* rawPtr = textureCoords->ptr();
const std::vector<double>* nncResultVals = m_nncData->connectionScalarResult(scalarResultIndex);
if (!nncResultVals)
{
textureCoords->setAll(cvf::Vec2f(0.0f, 1.0f));
return;
}
double cellScalarValue;
cvf::Vec2f texCoord;
@ -157,7 +164,7 @@ void RivNNCGeometryGenerator::textureCoordinates(cvf::Vec2fArray* textureCoords,
#pragma omp parallel for private(texCoord, cellScalarValue)
for (int tIdx = 0; tIdx < static_cast<int>(m_triangleIndexToNNCIndex->size()); tIdx++)
{
cellScalarValue = m_nncData->connections()[(*m_triangleIndexToNNCIndex)[tIdx]].m_transmissibility;
cellScalarValue = (*nncResultVals)[(*m_triangleIndexToNNCIndex)[tIdx]];
texCoord = mapper->mapToTextureCoord(cellScalarValue);
if (cellScalarValue == HUGE_VAL || cellScalarValue != cellScalarValue) // a != a is true for NAN's
{

View File

@ -44,7 +44,7 @@ public:
void setCellVisibility( const cvf::UByteArray* cellVisibilities, const RigGridBase * grid);
void textureCoordinates(cvf::Vec2fArray* textureCoords,
const cvf::ScalarMapper* mapper) const;
const cvf::ScalarMapper* mapper, size_t scalarResultIndex) const;
// Mapping between cells and geometry
cvf::ref<cvf::Array<size_t> > triangleToNNCIndex() const;

View File

@ -545,6 +545,10 @@ QString RimCase::relocateFile(const QString& orgFileName, const QString& orgNew
//--------------------------------------------------------------------------------------------------
bool RimCase::openReserviorCase()
{
// If read already, return
if (this->reservoirData() != NULL) return true;
if (!openEclipseGridFile())
{
return false;
@ -552,7 +556,14 @@ bool RimCase::openReserviorCase()
{
RimReservoirCellResultsStorage* results = this->results(RifReaderInterface::MATRIX_RESULTS);
if (results->cellResults()) results->cellResults()->createPlaceholderResultEntries();
if (results->cellResults())
{
results->cellResults()->createPlaceholderResultEntries();
// After the placeholder result for combined transmissibility is created,
// make sure the nnc transmissibilities can be addressed by this scalarResultIndex as well
size_t combinedTransResIdx = results->cellResults()->findScalarResultIndex(RimDefines::STATIC_NATIVE, RimDefines::combinedTransmissibilityResultName());
reservoirData()->mainGrid()->nncData()->setCombTransmisibilityScalarResultIndex(combinedTransResIdx);
}
}
{

View File

@ -342,6 +342,7 @@ size_t RimReservoirCellResultsStorage::findOrLoadScalarResult(RimDefines::Result
computeRiTransComponent(RimDefines::riTransXResultName());
computeRiTransComponent(RimDefines::riTransYResultName());
computeRiTransComponent(RimDefines::riTransZResultName());
computeNncCombRiTrans();
}
if ( resultName == RimDefines::riTransXResultName()
@ -438,8 +439,7 @@ void RimReservoirCellResultsStorage::computeSOILForTimeStep(size_t timeStepIndex
soilTimeStepCount = qMax(soilTimeStepCount, sgasTimeStepCount);
}
}
// Make sure memory is allocated for the new SOIL results
size_t soilResultScalarIndex = m_cellResults->findScalarResultIndex(RimDefines::DYNAMIC_NATIVE, "SOIL");
@ -453,7 +453,6 @@ void RimReservoirCellResultsStorage::computeSOILForTimeStep(size_t timeStepIndex
m_cellResults->cellScalarResults(soilResultScalarIndex, timeStepIndex).resize(soilResultValueCount);
std::vector<double>* swatForTimeStep = NULL;
std::vector<double>* sgasForTimeStep = NULL;
@ -735,8 +734,8 @@ void RimReservoirCellResultsStorage::computeRiTransComponent(const QString& riTr
// Get the needed result indices we depend on
size_t permResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, permCompName);
size_t ntgResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "NTG");
size_t ntgResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "NTG");
// Get the result index of the output
size_t riTransResultIdx = m_cellResults->findScalarResultIndex(RimDefines::STATIC_NATIVE, riTransComponentResultName);
@ -751,8 +750,8 @@ void RimReservoirCellResultsStorage::computeRiTransComponent(const QString& riTr
// Get all the actual result values
std::vector<double> & permResults = m_cellResults->cellScalarResults(permResultIdx)[0];
std::vector<double> & ntgResults = m_cellResults->cellScalarResults(ntgResultIdx)[0];
std::vector<double> & permResults = m_cellResults->cellScalarResults(permResultIdx)[0];
std::vector<double> & ntgResults = m_cellResults->cellScalarResults(ntgResultIdx)[0];
std::vector<double> & riTransResults = m_cellResults->cellScalarResults(riTransResultIdx)[0];
// Set up output container to correct number of results
@ -761,15 +760,15 @@ void RimReservoirCellResultsStorage::computeRiTransComponent(const QString& riTr
// Prepare how to index the result values:
bool isPermUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(permResultIdx);
bool isNtgUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(ntgResultIdx);
bool isPermUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(permResultIdx);
bool isNtgUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(ntgResultIdx);
bool isTransUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(riTransResultIdx);
// Set up result index function pointers
ResultIndexFunction riTranIdxFunc = isTransUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction permIdxFunc = isPermUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction ntgIdxFunc = isNtgUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction riTranIdxFunc = isTransUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction permIdxFunc = isPermUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction ntgIdxFunc = isNtgUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
const RigActiveCellInfo* activeCellInfo = m_cellResults->activeCellInfo();
@ -800,20 +799,20 @@ void RimReservoirCellResultsStorage::computeRiTransComponent(const QString& riTr
// Do nothing if neighbor cell has no results
size_t neighborCellPermResIdx = (*permIdxFunc)(activeCellInfo, neighborResvCellIdx);
if (neighborCellPermResIdx == cvf::UNDEFINED_SIZE_T) continue;
// Connection geometry
const RigFault* fault = grid->mainGrid()->findFaultFromCellIndexAndCellFace(nativeResvCellIndex, faceId);
bool isOnFault = fault;
cvf::Vec3d faceAreaVec;
cvf::Vec3d faceCenter;
if (isOnFault)
{
calculateConnectionGeometry( nativeCell, neighborCell, nodes,faceId,
&faceCenter, &faceAreaVec);
calculateConnectionGeometry(nativeCell, neighborCell, nodes, faceId,
&faceCenter, &faceAreaVec);
}
else
{
@ -863,7 +862,169 @@ void RimReservoirCellResultsStorage::computeRiTransComponent(const QString& riTr
}
}
}
#if 1
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimReservoirCellResultsStorage::computeNncCombRiTrans()
{
if (!m_cellResults) return;
size_t riCombTransScalarResultIndex = m_cellResults->findScalarResultIndex(RimDefines::STATIC_NATIVE, RimDefines::combinedRiTransResultName());
if (m_ownerMainGrid->nncData()->connectionScalarResult(riCombTransScalarResultIndex)) return;
// Todo: Get the correct one from Unit set read by ERT
double cdarchy = 0.008527; // (ECLIPSE 100) (METRIC)
// Get the needed result indices we depend on
size_t permXResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMX");
size_t permYResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMY");
size_t permZResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMZ");
size_t ntgResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "NTG");
// Get the result count, to handle that one of them might be globally defined
size_t permxResultValueCount = m_cellResults->cellScalarResults(permXResultIdx)[0].size();
size_t ntgResultValueCount = m_cellResults->cellScalarResults(ntgResultIdx)[0].size();
size_t resultValueCount = CVF_MIN(permxResultValueCount, ntgResultValueCount);
// Get all the actual result values
std::vector<double> & permXResults = m_cellResults->cellScalarResults(permXResultIdx)[0];
std::vector<double> & permYResults = m_cellResults->cellScalarResults(permYResultIdx)[0];
std::vector<double> & permZResults = m_cellResults->cellScalarResults(permZResultIdx)[0];
std::vector<double> & ntgResults = m_cellResults->cellScalarResults(ntgResultIdx)[0];
std::vector<double> & riCombTransResults = m_ownerMainGrid->nncData()->makeConnectionScalarResult(riCombTransScalarResultIndex);
// Prepare how to index the result values:
bool isPermXUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(permXResultIdx);
bool isPermYUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(permYResultIdx);
bool isPermZUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(permZResultIdx);
bool isNtgUsingResIdx = m_cellResults->isUsingGlobalActiveIndex(ntgResultIdx);
// Set up result index function pointers
ResultIndexFunction permXIdxFunc = isPermXUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction permYIdxFunc = isPermYUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction permZIdxFunc = isPermZUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
ResultIndexFunction ntgIdxFunc = isNtgUsingResIdx ? &reservoirActiveCellIndex : &directReservoirCellIndex;
const RigActiveCellInfo* activeCellInfo = m_cellResults->activeCellInfo();
const std::vector<cvf::Vec3d>& nodes = m_ownerMainGrid->nodes();
bool isFaceNormalsOutwards = m_ownerMainGrid->faceNormalsIsOutwards();
// NNC calculation
std::vector<RigConnection>& nncConnections = m_ownerMainGrid->nncData()->connections();
for (size_t connIdx = 0; connIdx < nncConnections.size(); connIdx++)
{
size_t nativeResvCellIndex = nncConnections[connIdx].m_c1GlobIdx;
size_t neighborResvCellIdx = nncConnections[connIdx].m_c2GlobIdx;
cvf::StructGridInterface::FaceType faceId = nncConnections[connIdx].m_c1Face;
ResultIndexFunction permIdxFunc = NULL;
std::vector<double> * permResults;
switch (faceId)
{
case cvf::StructGridInterface::POS_I:
case cvf::StructGridInterface::NEG_I:
permIdxFunc = permXIdxFunc;
permResults = &permXResults;
break;
case cvf::StructGridInterface::POS_J:
case cvf::StructGridInterface::NEG_J:
permIdxFunc = permYIdxFunc;
permResults = &permYResults;
break;
case cvf::StructGridInterface::POS_K:
case cvf::StructGridInterface::NEG_K:
permIdxFunc = permZIdxFunc;
permResults = &permZResults;
break;
}
// Do nothing if we are only dealing with active cells, and this cell is not active:
size_t nativeCellPermResIdx = (*permIdxFunc)(activeCellInfo, nativeResvCellIndex);
if (nativeCellPermResIdx == cvf::UNDEFINED_SIZE_T) continue;
// Do nothing if neighbor cell has no results
size_t neighborCellPermResIdx = (*permIdxFunc)(activeCellInfo, neighborResvCellIdx);
if (neighborCellPermResIdx == cvf::UNDEFINED_SIZE_T) continue;
const RigCell& nativeCell = m_ownerMainGrid->cells()[nativeResvCellIndex];
const RigCell& neighborCell = m_ownerMainGrid->cells()[neighborResvCellIdx];
// Connection geometry
cvf::Vec3d faceAreaVec = cvf::Vec3d::ZERO;;
cvf::Vec3d faceCenter = cvf::Vec3d::ZERO;;
// Polygon center
const std::vector<cvf::Vec3d>& realPolygon = nncConnections[connIdx].m_polygon;
for (size_t pIdx = 0; pIdx < realPolygon.size(); ++pIdx)
{
faceCenter += realPolygon[pIdx];
}
faceCenter *= 1.0 / realPolygon.size();
// Polygon area vector
faceAreaVec = cvf::GeometryTools::polygonAreaNormal3D(realPolygon);
if (!isFaceNormalsOutwards) faceAreaVec = -faceAreaVec;
double halfCellTrans = 0;
double neighborHalfCellTrans = 0;
// Native cell half cell transm
{
cvf::Vec3d centerToFace = faceCenter - nativeCell.center();
double perm = (*permResults)[nativeCellPermResIdx];
double ntg = 1.0;
if (faceId != cvf::StructGridInterface::POS_K)
{
size_t ntgResIdx = (*ntgIdxFunc)(activeCellInfo, nativeResvCellIndex);
ntg = ntgResults[ntgResIdx];
}
halfCellTrans = halfCellTransmissibility(perm, ntg, centerToFace, faceAreaVec);
}
// Neighbor cell half cell transm
{
cvf::Vec3d centerToFace = faceCenter - neighborCell.center();
double perm = (*permResults)[neighborCellPermResIdx];
double ntg = 1.0;
if (faceId != cvf::StructGridInterface::POS_K)
{
size_t ntgResIdx = (*ntgIdxFunc)(activeCellInfo, neighborResvCellIdx);
ntg = ntgResults[ntgResIdx];
}
neighborHalfCellTrans = halfCellTransmissibility(perm, ntg, centerToFace, -faceAreaVec);
}
double newtranTemp = newtran(cdarchy, 1.0, halfCellTrans, neighborHalfCellTrans);
riCombTransResults[connIdx] = newtranTemp;
}
#endif
}
#endif

View File

@ -66,10 +66,10 @@ protected:
private:
void computeSOILForTimeStep(size_t timeStepIndex);
void computeRiTransComponent(const QString& riTransComponentResultName);
void computeNncCombRiTrans();
QString getValidCacheFileName();
QString getCacheDirectoryPath();
// Fields
caf::PdmField<QString> m_resultCacheFileName;
caf::PdmPointersField<RimReservoirCellResultsStorageEntryInfo*>

View File

@ -1817,9 +1817,16 @@ void RimReservoirView::appendNNCResultInfo(size_t nncIndex, QString* resultInfo)
const RigConnection& conn = nncData->connections()[nncIndex];
cvf::StructGridInterface::FaceEnum face(conn.m_c1Face);
// Print result value for the NNC
size_t scalarResultIdx = this->cellResult()->gridScalarIndex();
const std::vector<double>* nncValues = nncData->connectionScalarResult(scalarResultIdx);
if (nncValues)
{
resultInfo->append(QString("NNC Value : %1\n").arg((*nncValues)[nncIndex]));
}
QString faultName;
resultInfo->append(QString("NNC Transmissibility : %1\n").arg(conn.m_transmissibility));
{
CVF_ASSERT(conn.m_c1GlobIdx < grid->cells().size());
const RigCell& cell = grid->cells()[conn.m_c1GlobIdx];

View File

@ -36,8 +36,7 @@ public:
RigConnection( )
: m_c1GlobIdx(cvf::UNDEFINED_SIZE_T),
m_c1Face(cvf::StructGridInterface::NO_FACE),
m_c2GlobIdx(cvf::UNDEFINED_SIZE_T),
m_transmissibility(0.0)
m_c2GlobIdx(cvf::UNDEFINED_SIZE_T)
{}
bool hasCommonArea() const
@ -49,8 +48,6 @@ public:
cvf::StructGridInterface::FaceType m_c1Face;
size_t m_c2GlobIdx;
double m_transmissibility;
std::vector<cvf::Vec3d> m_polygon;
};
@ -66,6 +63,33 @@ public:
std::vector<RigConnection>& connections() { return m_connections; }
const std::vector<RigConnection>& connections() const { return m_connections; };
std::vector<double>& makeConnectionScalarResult(size_t scalarResultIndex)
{
std::vector<double>& results = m_connectionResults[scalarResultIndex];
results.resize(m_connections.size(), HUGE_VAL);
return results;
}
const std::vector<double>* connectionScalarResult(size_t scalarResultIndex) const
{
std::map<size_t, std::vector<double> >::const_iterator it = m_connectionResults.find(scalarResultIndex);
if (it != m_connectionResults.end())
return &(it->second);
else
return NULL;
}
void setCombTransmisibilityScalarResultIndex(size_t scalarResultIndex)
{
std::map<size_t, std::vector<double> >::iterator it = m_connectionResults.find(cvf::UNDEFINED_SIZE_T);
CVF_ASSERT(it != m_connectionResults.end());
std::vector<double>& emptyData = m_connectionResults[scalarResultIndex];
std::vector<double>& realData = m_connectionResults[cvf::UNDEFINED_SIZE_T];
emptyData.swap(realData);
m_connectionResults.erase(cvf::UNDEFINED_SIZE_T);
}
private: // This section is possibly not needed
//const std::vector<size_t>& findConnectionIndices(size_t reservoirCellIndex, cvf::StructGridInterface::FaceType face) const;
//typedef std::map<size_t, caf::FixedArray<std::vector<size_t>, 7 > > ConnectionSearchMap;
@ -73,4 +97,5 @@ private: // This section is possibly not needed
private:
std::vector<RigConnection> m_connections;
std::map<size_t, std::vector<double> > m_connectionResults; ///< scalarResultIndex to value array map
};

View File

@ -541,6 +541,13 @@ void RigReservoirBuilderMock::addFaults(RigCaseData* eclipseCase)
addNnc(grid, i1, j1, k1, i2, j2, k2, nncConnections);
}
std::vector<double>& tranVals = grid->nncData()->makeConnectionScalarResult(cvf::UNDEFINED_SIZE_T);
for (size_t cIdx = 0; cIdx < tranVals.size(); ++cIdx)
{
tranVals[cIdx] = 0.2;
}
}
//--------------------------------------------------------------------------------------------------
@ -562,7 +569,6 @@ void RigReservoirBuilderMock::addNnc(RigMainGrid* grid, size_t i1, size_t j1, si
RigConnection conn;
conn.m_c1GlobIdx = c1GlobalIndex;
conn.m_c2GlobIdx = c2GlobalIndex;
conn.m_transmissibility = 0.2;
nncConnections.push_back(conn);
}