From 57a049d39f92f23a4cee6bdd2180fa912f130633 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Fri, 22 Aug 2014 08:07:00 +0200 Subject: [PATCH] riTRANSXYZ: Moving the calculation to ReservoirCellResultsStorage Create stored native riTRANX,Y,Z results instead of doing everyting in one go. This makes the statistics calculation work automatically instead of beeing custoized and require a compplete calculation of its own. --- .../RimReservoirCellResultsStorage.cpp | 241 ++++++++++++++++++ .../RimReservoirCellResultsStorage.h | 2 +- .../RigCaseCellResultsData.h | 5 + .../ReservoirDataModel/RigCaseData.cpp | 5 + .../RigCombRiTransResultAccessor.cpp | 27 +- .../RigCombRiTransResultAccessor.h | 4 +- 6 files changed, 269 insertions(+), 15 deletions(-) diff --git a/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.cpp b/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.cpp index c4f2e8a7a5..5baa428661 100644 --- a/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.cpp +++ b/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.cpp @@ -19,7 +19,9 @@ #include "RimReservoirCellResultsStorage.h" #include "RigCaseCellResultsData.h" +#include "RigActiveCellInfo.h" #include "RigMainGrid.h" +#include "RigCell.h" #include "RimTools.h" #include "cafProgressInfo.h" @@ -29,6 +31,7 @@ #include #include #include +#include "cvfGeometryTools.h" CAF_PDM_SOURCE_INIT(RimReservoirCellResultsStorage, "ReservoirCellResultStorage"); @@ -564,6 +567,242 @@ void RimReservoirCellResultsStorage::computeDepthRelatedResults() } } +#if 1 + +void calculateConnectionGeometry(const RigCell& c1, const RigCell& c2 , const std::vector& nodes, + cvf::StructGridInterface::FaceType faceId, + cvf::Vec3d* faceCenter, cvf::Vec3d* faceAreaVec) +{ + CVF_TIGHT_ASSERT(faceCenter && faceAreaVec); + + *faceCenter = cvf::Vec3d::ZERO; + *faceAreaVec = cvf::Vec3d::ZERO; + + std::vector polygon; + std::vector intersections; + caf::SizeTArray4 face1; + caf::SizeTArray4 face2; + c1.faceIndices(faceId, &face1); + c2.faceIndices(cvf::StructGridInterface::oppositeFace(faceId), &face2); + + bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( + &polygon, + &intersections, + (cvf::EdgeIntersectStorage*)NULL, + cvf::wrapArrayConst(&nodes), + face1.data(), + face2.data(), + 1e-6); + + + if (foundOverlap) + { + std::vector realPolygon; + + for (size_t pIdx = 0; pIdx < polygon.size(); ++pIdx) + { + if (polygon[pIdx] < nodes.size()) + realPolygon.push_back(nodes[polygon[pIdx]]); + else + realPolygon.push_back(intersections[polygon[pIdx] - nodes.size()]); + } + + // Polygon center + 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); + + } + +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double halfCellTransmissibility(double perm, double ntg, const cvf::Vec3d& centerToFace, const cvf::Vec3d& faceAreaVec) +{ + return perm*ntg*(faceAreaVec*centerToFace) / (centerToFace*centerToFace); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double newtran(double cdarchy, double mult, double halfCellTrans, double neighborHalfCellTrans) +{ + return cdarchy * mult / ((1 / halfCellTrans) + (1 / neighborHalfCellTrans)); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +typedef size_t (*ResultIndexFunction)(const RigActiveCellInfo* activeCellinfo, size_t reservoirCellIndex); + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- + +size_t directReservoirCellIndex(const RigActiveCellInfo* activeCellinfo, size_t reservoirCellIndex) +{ + return reservoirCellIndex; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- + +size_t reservoirActiveCellIndex(const RigActiveCellInfo* activeCellinfo, size_t reservoirCellIndex) +{ + return activeCellinfo->cellResultIndex(reservoirCellIndex); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimReservoirCellResultsStorage::computeRiTransX() +{ + if (!m_cellResults) return; + + cvf::StructGridInterface::FaceType faceId = cvf::StructGridInterface::POS_I; + QString riTRANName = "riTRANX"; + QString permName = "PERMX"; + + double cdarchy = 0.008527; // (ECLIPSE 100) (METRIC) + + // Get the needed result indices we depend on + + size_t permResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMX"); + 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(permResultIdx)[0].size(); + size_t ntgResultValueCount = m_cellResults->cellScalarResults(ntgResultIdx)[0].size(); + + size_t resultValueCount = CVF_MIN(permxResultValueCount, ntgResultValueCount); + + // Check if we need to create a placeholder result, or if it exists already + + size_t riTransResultIdx = findOrLoadScalarResult(RimDefines::STATIC_NATIVE, riTRANName); + + if (riTransResultIdx == cvf::UNDEFINED_SIZE_T) + { + riTransResultIdx = m_cellResults->addStaticScalarResult(RimDefines::STATIC_NATIVE, riTRANName, false, resultValueCount); + } + + // Prepare how to index the result values: + + 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; + + // Get all the actual result values + + std::vector & riTransResults = m_cellResults->cellScalarResults(riTransResultIdx)[0]; + std::vector & permResults = m_cellResults->cellScalarResults(permResultIdx)[0]; + std::vector & ntgResults = m_cellResults->cellScalarResults(ntgResultIdx)[0]; + + const RigActiveCellInfo* activeCellInfo = m_cellResults->activeCellInfo(); + const std::vector& nodes = m_ownerMainGrid->nodes(); + + for (size_t nativeResvCellIndex = 0; nativeResvCellIndex < m_ownerMainGrid->cells().size(); nativeResvCellIndex++) + { + // Do nothing if we are only dealing with active cells, and this cell is not active: + if ((*riTranIdxFunc)(activeCellInfo, nativeResvCellIndex) == cvf::UNDEFINED_SIZE_T) continue; + + const RigCell& nativeCell = m_ownerMainGrid->cells()[nativeResvCellIndex]; + RigGridBase* grid = nativeCell.hostGrid(); + + size_t gridLocalNativeCellIndex = nativeCell.gridLocalCellIndex(); + + size_t i, j, k, gridLocalNeighborCellIdx; + + grid->ijkFromCellIndex(gridLocalNativeCellIndex, &i, &j, &k); + + if (grid->cellIJKNeighbor(i, j, k, faceId, &gridLocalNeighborCellIdx)) + { + size_t neighborResvCellIdx = grid->reservoirCellIndex(gridLocalNeighborCellIdx); + const RigCell& neighborCell = m_ownerMainGrid->cells()[neighborResvCellIdx]; + + // 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); + } + else + { + faceCenter = nativeCell.faceCenter(faceId); + faceAreaVec = nativeCell.faceNormalWithAreaLenght(faceId); + } + + double halfCellTrans = 0; + double neighborHalfCellTrans = 0; + + // Native cell half cell transm + { + cvf::Vec3d centerToFace = faceCenter - nativeCell.center(); + + size_t permResIdx = (*permIdxFunc)(activeCellInfo, nativeResvCellIndex); + double perm = permResults[permResIdx]; + + 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(); + + size_t permResIdx = (*permIdxFunc)(activeCellInfo, neighborResvCellIdx); + double perm = permResults[permResIdx]; + + 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); + } + + size_t tranResIdx = (*riTranIdxFunc)(activeCellInfo, nativeResvCellIndex); + riTransResults[tranResIdx] = newtran(cdarchy, 1.0, halfCellTrans, neighborHalfCellTrans);; + } + + } + +} + +#endif //-------------------------------------------------------------------------------------------------- /// @@ -594,6 +833,8 @@ size_t RimReservoirCellResultsStorage::findOrLoadScalarResult(const QString& res return scalarResultIndex; } + + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.h b/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.h index 55713ce290..78631bfe41 100644 --- a/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.h +++ b/ApplicationCode/ProjectDataModel/RimReservoirCellResultsStorage.h @@ -69,7 +69,7 @@ private: QString getValidCacheFileName(); QString getCacheDirectoryPath(); - + void computeRiTransX(); // Fields caf::PdmField m_resultCacheFileName; caf::PdmPointersField diff --git a/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h b/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h index b236bc4b1b..30cf5c3ac6 100644 --- a/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h +++ b/ApplicationCode/ReservoirDataModel/RigCaseCellResultsData.h @@ -30,6 +30,7 @@ class RifReaderInterface; class RigMainGrid; class RigStatisticsDataCache; +class RigActiveCellInfo; //================================================================================================== /// Class containing the results for the complete number of active cells. Both main grid and LGR's @@ -40,6 +41,9 @@ public: RigCaseCellResultsData(RigMainGrid* ownerGrid); void setMainGrid(RigMainGrid* ownerGrid); + void setActiveCellInfo(RigActiveCellInfo* activeCellInfo) { m_activeCellInfo = activeCellInfo;} + RigActiveCellInfo* activeCellInfo() { return m_activeCellInfo;} + const RigActiveCellInfo* activeCellInfo() const { return m_activeCellInfo;} // Max and min values of the results void recalculateStatistics(size_t scalarResultIndex); @@ -123,5 +127,6 @@ private: std::vector m_resultInfos; RigMainGrid* m_ownerMainGrid; + RigActiveCellInfo* m_activeCellInfo; }; diff --git a/ApplicationCode/ReservoirDataModel/RigCaseData.cpp b/ApplicationCode/ReservoirDataModel/RigCaseData.cpp index 2d7b204c48..3dd84a69b8 100644 --- a/ApplicationCode/ReservoirDataModel/RigCaseData.cpp +++ b/ApplicationCode/ReservoirDataModel/RigCaseData.cpp @@ -33,6 +33,9 @@ RigCaseData::RigCaseData() m_activeCellInfo = new RigActiveCellInfo; m_fractureActiveCellInfo = new RigActiveCellInfo; + + m_matrixModelResults->setActiveCellInfo(m_activeCellInfo.p()); + m_fractureModelResults->setActiveCellInfo(m_fractureActiveCellInfo.p()); } //-------------------------------------------------------------------------------------------------- @@ -397,10 +400,12 @@ void RigCaseData::setActiveCellInfo(RifReaderInterface::PorosityModelResultType if (porosityModel == RifReaderInterface::MATRIX_RESULTS) { m_activeCellInfo = activeCellInfo; + m_matrixModelResults->setActiveCellInfo(m_activeCellInfo.p()); } else { m_fractureActiveCellInfo = activeCellInfo; + m_fractureModelResults->setActiveCellInfo(m_fractureActiveCellInfo.p()); } } diff --git a/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.cpp b/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.cpp index 78b72b8d34..5f95b1eb1f 100644 --- a/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.cpp +++ b/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.cpp @@ -112,7 +112,7 @@ double RigCombRiTransResultAccessor::getNtgValue(size_t gridLocalCellIndex, cvf: //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double halfCellTransmissibility( double perm , double ntg, const cvf::Vec3d& centerToFace, const cvf::Vec3d& faceAreaVec) +double RigCombRiTransResultAccessor::halfCellTransmissibility(double perm, double ntg, const cvf::Vec3d& centerToFace, const cvf::Vec3d& faceAreaVec) { return perm*ntg*(faceAreaVec*centerToFace) / (centerToFace*centerToFace); } @@ -120,11 +120,12 @@ double halfCellTransmissibility( double perm , double ntg, const cvf::Vec3d& cen //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double newtran(double cdarchy, double mult, double halfCellTrans, double neighborHalfCellTrans) +double RigCombRiTransResultAccessor::newtran(double cdarchy, double mult, double halfCellTrans, double neighborHalfCellTrans) { return cdarchy * mult / ( ( 1 / halfCellTrans) + (1 / neighborHalfCellTrans) ); } + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -139,7 +140,7 @@ double RigCombRiTransResultAccessor::calculateHalfCellTrans(size_t gridLocalCell if (isFaultFace) { - calculateConnectionGeometry( gridLocalCellIndex, neighborGridCellIdx, faceId, &faceCenter, &faceAreaVec); + calculateConnectionGeometry( m_grid, gridLocalCellIndex, neighborGridCellIdx, faceId, &faceCenter, &faceAreaVec); } else { @@ -157,18 +158,18 @@ double RigCombRiTransResultAccessor::calculateHalfCellTrans(size_t gridLocalCell return halfCellTransmissibility(perm, ntg, centerToFace, faceAreaVec); } -void RigCombRiTransResultAccessor::calculateConnectionGeometry( size_t gridLocalCellIndex, size_t neighborGridCellIdx, +void RigCombRiTransResultAccessor::calculateConnectionGeometry( const RigGridBase* grid, size_t gridLocalCellIndex, size_t neighborGridCellIdx, cvf::StructGridInterface::FaceType faceId, - cvf::Vec3d* faceCenter, cvf::Vec3d* faceAreaVec) const + cvf::Vec3d* faceCenter, cvf::Vec3d* faceAreaVec) { CVF_TIGHT_ASSERT(faceCenter && faceAreaVec); *faceCenter = cvf::Vec3d::ZERO; *faceAreaVec = cvf::Vec3d::ZERO; - const RigMainGrid* mainGrid = m_grid->mainGrid(); + const RigMainGrid* mainGrid = grid->mainGrid(); - const RigCell& c1 = m_grid->cell(gridLocalCellIndex); - const RigCell& c2 = m_grid->cell(neighborGridCellIdx); + const RigCell& c1 = grid->cell(gridLocalCellIndex); + const RigCell& c2 = grid->cell(neighborGridCellIdx); std::vector polygon; std::vector intersections; @@ -193,10 +194,10 @@ void RigCombRiTransResultAccessor::calculateConnectionGeometry( size_t gridLocal for (size_t pIdx = 0; pIdx < polygon.size(); ++pIdx) { - if (polygon[pIdx] < m_grid->mainGrid()->nodes().size()) - realPolygon.push_back(m_grid->mainGrid()->nodes()[polygon[pIdx]]); + if (polygon[pIdx] < grid->mainGrid()->nodes().size()) + realPolygon.push_back(grid->mainGrid()->nodes()[polygon[pIdx]]); else - realPolygon.push_back(intersections[polygon[pIdx] - m_grid->mainGrid()->nodes().size()]); + realPolygon.push_back(intersections[polygon[pIdx] - grid->mainGrid()->nodes().size()]); } // Polygon center @@ -216,8 +217,8 @@ void RigCombRiTransResultAccessor::calculateConnectionGeometry( size_t gridLocal } //-------------------------------------------------------------------------------------------------- /// -/// Neighbour cell transmisibilities only is calculated here -/// Not NNC transmisibilities. Thart has to be done separatly elsewhere +/// Neighbor cell transmisibilities only is calculated here +/// Not NNC transmisibilities. That has to be done separately elsewhere /// /// Todo: What about Grid to Grid connections ? /// Todo: needs optimization. Things are done several times. Caching of the results should be considered. etc. diff --git a/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.h b/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.h index 7576272976..9ffe227112 100644 --- a/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.h +++ b/ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.h @@ -53,8 +53,10 @@ private: double calculateHalfCellTrans( size_t gridLocalCellIndex, size_t neighborGridCellIdx, cvf::StructGridInterface::FaceType faceId, bool isFaultFace) const; double getNtgValue( size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId ) const; - void calculateConnectionGeometry( size_t gridLocalCellIndex, size_t neighborGridCellIdx, cvf::StructGridInterface::FaceType faceId, cvf::Vec3d* centerToFace, cvf::Vec3d* faceAreaVec) const; + static void calculateConnectionGeometry( const RigGridBase* grid, size_t gridLocalCellIndex, size_t neighborGridCellIdx, cvf::StructGridInterface::FaceType faceId, cvf::Vec3d* centerToFace, cvf::Vec3d* faceAreaVec); + static double halfCellTransmissibility(double perm, double ntg, const cvf::Vec3d& centerToFace, const cvf::Vec3d& faceAreaVec); + static double newtran(double cdarchy, double mult, double halfCellTrans, double neighborHalfCellTrans); cvf::ref m_xPermAccessor; cvf::ref m_yPermAccessor;