diff --git a/ApplicationCode/ProjectDataModel/RimStimPlanCell.h b/ApplicationCode/ProjectDataModel/RimStimPlanCell.h index a9ee97f628..92ec869e19 100644 --- a/ApplicationCode/ProjectDataModel/RimStimPlanCell.h +++ b/ApplicationCode/ProjectDataModel/RimStimPlanCell.h @@ -42,8 +42,8 @@ public: std::vector getPolygon() { return m_polygon; } double getValue() { return m_value; } - size_t i() { return m_i; } - size_t j() { return m_j; } + size_t getI() { return m_i; } + size_t getJ() { return m_j; } private: diff --git a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp index 989c35c948..1e4cf4e721 100644 --- a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp +++ b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp @@ -804,6 +804,52 @@ std::vector RimStimPlanFractureTemplate::getStimPlanCells(cons return stimPlanCells; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimStimPlanFractureTemplate::getStimPlanRowPolygon(size_t i) //Row with constant depth +{ + std::vector rowPolygon; + + std::vector depthCoordsAtNodes = adjustedDepthCoordsAroundWellPathPosition(); + std::vector xCoordsAtNodes = getNegAndPosXcoords(); + + std::vector xCoords; + for (int i = 0; i < xCoordsAtNodes.size() - 1; i++) xCoords.push_back((xCoordsAtNodes[i] + xCoordsAtNodes[i + 1]) / 2); + std::vector depthCoords; + for (int i = 0; i < depthCoordsAtNodes.size() - 1; i++) depthCoords.push_back((depthCoordsAtNodes[i] + depthCoordsAtNodes[i + 1]) / 2); + + rowPolygon.push_back(cvf::Vec3d(static_cast(xCoords[0]), static_cast(depthCoords[i]), 0.0)); + rowPolygon.push_back(cvf::Vec3d(static_cast(xCoords.back()), static_cast(depthCoords[i]), 0.0)); + rowPolygon.push_back(cvf::Vec3d(static_cast(xCoords.back()), static_cast(depthCoords[i+1]), 0.0)); + rowPolygon.push_back(cvf::Vec3d(static_cast(xCoords[0]), static_cast(depthCoords[i+1]), 0.0)); + + return rowPolygon; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimStimPlanFractureTemplate::getStimPlanColPolygon(size_t j) +{ + std::vector colPolygon; + + std::vector depthCoordsAtNodes = adjustedDepthCoordsAroundWellPathPosition(); + std::vector xCoordsAtNodes = getNegAndPosXcoords(); + + std::vector xCoords; + for (int i = 0; i < xCoordsAtNodes.size() - 1; i++) xCoords.push_back((xCoordsAtNodes[i] + xCoordsAtNodes[i + 1]) / 2); + std::vector depthCoords; + for (int i = 0; i < depthCoordsAtNodes.size() - 1; i++) depthCoords.push_back((depthCoordsAtNodes[i] + depthCoordsAtNodes[i + 1]) / 2); + + colPolygon.push_back(cvf::Vec3d(static_cast(xCoords[j]), static_cast(depthCoords[0]), 0.0)); + colPolygon.push_back(cvf::Vec3d(static_cast(xCoords[j+1]), static_cast(depthCoords[0]), 0.0)); + colPolygon.push_back(cvf::Vec3d(static_cast(xCoords[j+1]), static_cast(depthCoords.back()), 0.0)); + colPolygon.push_back(cvf::Vec3d(static_cast(xCoords[j]), static_cast(depthCoords.back()), 0.0)); + + return colPolygon; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -926,6 +972,22 @@ void RimStimPlanFractureTemplate::sortPolygon(std::vector &polygon) } } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RimStimPlanFractureTemplate::stimPlanGridNumberOfColums() +{ + return getNegAndPosXcoords().size(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RimStimPlanFractureTemplate::stimPlanGridNumberOfRows() +{ + return adjustedDepthCoordsAroundWellPathPosition().size(); +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h index 7adea58dc5..79e4942f25 100644 --- a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h +++ b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h @@ -65,6 +65,10 @@ public: std::vector fracturePolygon(caf::AppEnum< RimDefines::UnitSystem > fractureUnit); void sortPolygon(std::vector &polygon); + size_t stimPlanGridNumberOfRows(); + size_t stimPlanGridNumberOfColums(); + + std::vector getNegAndPosXcoords(); std::vector adjustedDepthCoordsAroundWellPathPosition(); std::vector getStimPlanTimeValues(); @@ -73,6 +77,9 @@ public: void getStimPlanDataAsPolygonsAndValues(std::vector > &cellsAsPolygons, std::vector ¶meterValue, const QString& resultName, const QString& unitName, size_t timeStepIndex); std::vector getStimPlanCells(const QString& resultName, const QString& unitName, size_t timeStepIndex); + std::vector getStimPlanRowPolygon(size_t i); + std::vector getStimPlanColPolygon(size_t j); + void loadDataAndUpdate(); void setDefaultsBasedOnXMLfile(); diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp index 6ce3a49c17..e8f9172677 100644 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp @@ -38,6 +38,7 @@ #include "RigMainGrid.h" #include "cvfMath.h" #include "RimDefines.h" +#include "RimStimPlanCell.h" //-------------------------------------------------------------------------------------------------- /// @@ -404,6 +405,134 @@ void RigFractureTransCalc::computeUpscaledPropertyFromStimPlanForEclipseCell(dou } } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigFractureTransCalc::HAflowAcrossLayersUpscale(QString resultName, QString resultUnit, size_t timeStepIndex, RimDefines::UnitSystem unitSystem, size_t eclipseCellIndex) +{ + + //TODO: A lot of common code with function for calculating transmissibility... + + RimStimPlanFractureTemplate* fracTemplateStimPlan; + if (dynamic_cast(m_fracture->attachedFractureDefinition())) + { + fracTemplateStimPlan = dynamic_cast(m_fracture->attachedFractureDefinition()); + } + else return cvf::UNDEFINED_DOUBLE; + + std::vector stimPlanCells = fracTemplateStimPlan->getStimPlanCells(resultName, resultUnit, timeStepIndex); + +// RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS; +// RimReservoirCellResultsStorage* gridCellResults = m_case->results(porosityModel); + + + cvf::Vec3d localX; + cvf::Vec3d localY; + cvf::Vec3d localZ; + std::vector > planeCellPolygons; + bool isPlanIntersected = planeCellIntersectionPolygons(eclipseCellIndex, planeCellPolygons, localX, localY, localZ); + if (!isPlanIntersected || planeCellPolygons.size() == 0) return cvf::UNDEFINED_DOUBLE; + + //Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon/stimPlan mesh already is located) + cvf::Mat4f invertedTransMatrix = m_fracture->transformMatrix().getInverted(); + for (std::vector & planeCellPolygon : planeCellPolygons) + { + for (cvf::Vec3d& v : planeCellPolygon) + { + v.transformPoint(static_cast(invertedTransMatrix)); + } + } + + //Vector for vertical - på tvers av lag... Gitt av orientering av stimPlan-grid... + + cvf::Vec3d directionAcrossLayers; + cvf::Vec3d directionAlongLayers; + + //TODO: Get these vectors properly... + directionAcrossLayers = cvf::Vec3d(0, -1, 0); + directionAlongLayers = cvf::Vec3d(1, 0, 0); + + directionAcrossLayers.transformVector(static_cast(invertedTransMatrix)); + directionAlongLayers.transformVector(static_cast(invertedTransMatrix)); + + std::vector fracPolygon = m_fracture->attachedFractureDefinition()->fracturePolygon(unitSystem); + std::vector > polygonsDescribingFractureInCell; + + std::vector upscaledConductivities; + + + //Harmonic weighted mean + for (std::vector planeCellPolygon : planeCellPolygons) + { + std::vector DcolSum; + std::vector lavgCol; + std::vector CondHarmCol; + + for (size_t j = 0; j < fracTemplateStimPlan->stimPlanGridNumberOfColums(); j++) + { + std::vector conductivitiesInStimPlanCells; + std::vector lengthsLiOfStimPlanCol; + std::vector heightsDioFStimPlanCells; + + std::vector stimPlanCellsCol = getColOfStimPlanCells(stimPlanCells, j); + for (RimStimPlanCell* stimPlanCell : stimPlanCellsCol) + { + if (stimPlanCell->getValue() > 1e-7) + { + std::vector >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell->getPolygon(), planeCellPolygon); + if (clippedStimPlanPolygons.size() > 0) + { + for (auto clippedStimPlanPolygon : clippedStimPlanPolygons) + { + conductivitiesInStimPlanCells.push_back(stimPlanCell->getValue()); + lengthsLiOfStimPlanCol.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAlongLayers, clippedStimPlanPolygon)); + heightsDioFStimPlanCells.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAcrossLayers, clippedStimPlanPolygon)); + } + } + } + } + //Regne ut average + double sumDiDivCondLi = 0.0; + double sumDi = 0.0; + double sumLiDi = 0.0; + for (int i = 0; i < conductivitiesInStimPlanCells.size(); i++) + { + sumDiDivCondLi += heightsDioFStimPlanCells[i] / (conductivitiesInStimPlanCells[i] * lengthsLiOfStimPlanCol[i]); + sumDi += heightsDioFStimPlanCells[i]; + sumLiDi += heightsDioFStimPlanCells[i] * lengthsLiOfStimPlanCol[i]; + } + + if (sumDiDivCondLi != 0) + { + DcolSum.push_back(sumDi); + double lAvgValue = sumLiDi / sumDi; + lavgCol.push_back(lAvgValue); + CondHarmCol.push_back(1 / (sumLiDi*sumDiDivCondLi)); + } + } + + //Do arithmetic upscaling based on harmonic upscaled values for coloums + double sumCondHLiDivDi = 0.0; + double sumLi = 0.0; + double sumDiLi = 0.0; + for (int i = 0; i < CondHarmCol.size(); i++) + { + sumLi += lavgCol[i]; + sumDiLi += DcolSum[i] * lavgCol[i]; + sumCondHLiDivDi += CondHarmCol[i] * lavgCol[i] / DcolSum[i]; + } + double Davg = sumDiLi / sumLi; + double condHA = (Davg / sumLi) * sumCondHLiDivDi; + + upscaledConductivities.push_back(condHA); + + } + + //TODO: Is this the right way of handling getting several values for each cell? + return arithmeticAverage(upscaledConductivities); + +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -448,6 +577,25 @@ double RigFractureTransCalc::areaWeightedArithmeticAverage(std::vector a } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigFractureTransCalc::arithmeticAverage(std::vector values) +{ + if (values.size() == 0) return cvf::UNDEFINED_DOUBLE; + + double sumValue = 0.0; + size_t numberOfValues = 0; + for (double value : values) + { + sumValue += value; + numberOfValues++; + } + + return sumValue / numberOfValues; +} + + //-------------------------------------------------------------------------------------------------- /// @@ -600,6 +748,42 @@ void RigFractureTransCalc::computeFlowIntoTransverseWell() } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigFractureTransCalc::getRowOfStimPlanCells(std::vector allStimPlanCells, size_t i) +{ + std::vector stimPlanCellRow; + + for (RimStimPlanCell* stimPlanCell : allStimPlanCells) + { + if (stimPlanCell->getI() == i) + { + stimPlanCellRow.push_back(stimPlanCell); + } + } + + return stimPlanCellRow; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigFractureTransCalc::getColOfStimPlanCells(std::vector allStimPlanCells, size_t j) +{ + std::vector stimPlanCellCol; + + for (RimStimPlanCell* stimPlanCell : allStimPlanCells) + { + if (stimPlanCell->getJ() == j) + { + stimPlanCellCol.push_back(stimPlanCell); + } + } + + return stimPlanCellCol; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h index 6d841170fe..5212dcda7b 100644 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h @@ -35,7 +35,7 @@ class RimFracture; class RimEclipseCase; - +class RimStimPlanCell; //================================================================================================== /// @@ -52,13 +52,22 @@ public: void computeUpscaledPropertyFromStimPlan(QString resultName, QString resultUnit, size_t timeStepIndex); void computeUpscaledPropertyFromStimPlanForEclipseCell(double &upscaledAritmStimPlanValue, double &upscaledHarmStimPlanValue, QString resultName, QString resultUnit, size_t timeStepIndex, caf::AppEnum< RimDefines::UnitSystem > unitSystem, size_t cellIndex); + double HAflowAcrossLayersUpscale(QString resultName, QString resultUnit, size_t timeStepIndex, RimDefines::UnitSystem unitSystem, size_t eclipseCellIndex); + + + static double areaWeightedHarmonicAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts); static double areaWeightedArithmeticAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts); - + static double arithmeticAverage(std::vector values); + void computeFlowInFracture(); void computeFlowIntoTransverseWell(); + static std::vector getRowOfStimPlanCells(std::vector allStimPlanCells, size_t i); + static std::vector getColOfStimPlanCells(std::vector allStimPlanCells, size_t j); + + private: RimEclipseCase* m_case; RimFracture* m_fracture;