From ffeb50b06a7f8874ce7d4e608097c4957953fd58 Mon Sep 17 00:00:00 2001 From: astridkbjorke Date: Fri, 17 Mar 2017 11:28:00 +0100 Subject: [PATCH] #1320 - pre-proto - Prototyping calculation of upscaled parameter values from stimPlan. Bug in calculation so output values are wrong. --- .../FileInterface/RifEclipseExportTools.cpp | 40 ++++- .../RivWellFracturePartMgr.h | 3 +- .../ProjectDataModel/RimFracture.cpp | 147 ++++++++++++++++++ .../ProjectDataModel/RimFracture.h | 2 + .../RimStimPlanFractureTemplate.cpp | 61 +++++++- .../RimStimPlanFractureTemplate.h | 17 +- .../ReservoirDataModel/RigFracture.h | 1 + 7 files changed, 258 insertions(+), 13 deletions(-) diff --git a/ApplicationCode/FileInterface/RifEclipseExportTools.cpp b/ApplicationCode/FileInterface/RifEclipseExportTools.cpp index 7edefbba71..32219bcaef 100644 --- a/ApplicationCode/FileInterface/RifEclipseExportTools.cpp +++ b/ApplicationCode/FileInterface/RifEclipseExportTools.cpp @@ -83,7 +83,6 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c } caf::ProgressInfo pi(fractures.size(), QString("Writing data to file %1").arg(fileName)); - RimEclipseWell* simWell = nullptr; RimWellPath* wellPath = nullptr; @@ -99,6 +98,45 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c if (caseUnit == RigEclipseCaseData::UNITS_FIELD) out << "-- Using field unit system" << "\n"; out << "\n"; + + for (RimFracture* fracture : fractures) //For testing upscaling... + { + fracture->computeUpscaledPropertyFromStimPlan(caseToApply); + std::vector fracDataVector = fracture->attachedRigFracture()->fractureData(); + + for (RigFractureData fracData : fracDataVector) + { + out << qSetFieldWidth(4); + out << "-- "; + + out << qSetFieldWidth(12); + wellPath, simWell = nullptr; + fracture->firstAncestorOrThisOfType(simWell); + if (simWell) out << simWell->name + " "; // 1. Well name + fracture->firstAncestorOrThisOfType(wellPath); + if (wellPath) out << wellPath->name + " "; // 1. Well name + + out << qSetFieldWidth(16); + out << fracture->name().left(15) + " "; + + + out << qSetFieldWidth(5); + size_t i, j, k; + mainGrid->ijkFromCellIndex(fracData.reservoirCellIndex, &i, &j, &k); + out << i + 1; // 2. I location grid block, adding 1 to go to eclipse 1-based grid definition + out << j + 1; // 3. J location grid block, adding 1 to go to eclipse 1-based grid definition + out << k + 1; // 4. K location of upper connecting grid block, adding 1 to go to eclipse 1-based grid definition + + out << qSetFieldWidth(10); + out << QString::number(fracData.upscaledStimPlanValue, 'f', 3); + + out << "\n"; + } + } + + + + printBackgroundDataHeaderLine(out); diff --git a/ApplicationCode/ModelVisualization/RivWellFracturePartMgr.h b/ApplicationCode/ModelVisualization/RivWellFracturePartMgr.h index 414aec2165..0c4f5f8981 100644 --- a/ApplicationCode/ModelVisualization/RivWellFracturePartMgr.h +++ b/ApplicationCode/ModelVisualization/RivWellFracturePartMgr.h @@ -54,6 +54,8 @@ public: void appendGeometryPartsToModel(cvf::ModelBasicList* model, caf::DisplayCoordTransform* displayCoordTransform); void clearGeometryCache(); + static std::vector mirrorDataAtSingleDepth(std::vector depthData); + private: void updatePartGeometry(caf::DisplayCoordTransform* displayCoordTransform); void updatePartGeometryTexture(caf::DisplayCoordTransform* displayCoordTransform); @@ -69,7 +71,6 @@ private: std::vector transfromToFractureDisplayCoords(std::vector polygon, cvf::Mat4f m, caf::DisplayCoordTransform* displayCoordTransform); bool stimPlanCellTouchesPolygon(const std::vector& polygon, double xMin, double xMax, double yMin, double yMax, float polygonXmin, float polygonXmax, float polygonYmin, float polygonYmax); - std::vector mirrorDataAtSingleDepth(std::vector depthData); static cvf::ref createGeo(const std::vector& triangleIndices, const std::vector& nodeCoords); diff --git a/ApplicationCode/ProjectDataModel/RimFracture.cpp b/ApplicationCode/ProjectDataModel/RimFracture.cpp index b0ef6c463c..deae8b2a58 100644 --- a/ApplicationCode/ProjectDataModel/RimFracture.cpp +++ b/ApplicationCode/ProjectDataModel/RimFracture.cpp @@ -497,6 +497,153 @@ void RimFracture::computeTransmissibility(RimEclipseCase* caseToApply) } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimFracture::computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply) +{ + if (!attachedFractureDefinition()) return; + + RimStimPlanFractureTemplate* fracTemplateStimPlan; + if (dynamic_cast(attachedFractureDefinition())) + { + fracTemplateStimPlan = dynamic_cast(attachedFractureDefinition()); + } + else return; + + //Get correct unit system: + RigEclipseCaseData::UnitsType caseUnit = caseToApply->eclipseCaseData()->unitsType(); + RimDefines::UnitSystem unitForExport; + + if (caseUnit == RigEclipseCaseData::UNITS_METRIC) + { + RiaLogging::debug(QString("Calculating transmissibilities for %1 in metric units").arg(name())); + unitForExport = RimDefines::UNITS_METRIC; + } + else if (caseUnit == RigEclipseCaseData::UNITS_FIELD) + { + RiaLogging::debug(QString("Calculating transmissibilities for %1 in field units").arg(name())); + unitForExport = RimDefines::UNITS_FIELD; + } + else + { + RiaLogging::error(QString("Unit system for case not supported for fracture export.")); + return; + } + + + //TODO: Get these more generally: + QString resultName = "CONDUCTIVITY"; + QString resultUnit = "md-m"; + size_t timeStepIndex = 0; + + + std::vector > stimPlanCellsAsPolygons; + std::vector stimPlanParameterValues; + fracTemplateStimPlan->getStimPlanDataAsPolygonsAndValues(stimPlanCellsAsPolygons, stimPlanParameterValues, resultName, resultUnit, timeStepIndex); + + + + //TODO: A lot of common code with function above... Can be cleaned up...? + std::vector fracCells = getPotentiallyFracturedCells(); + + + RigEclipseCaseData* eclipseCaseData = caseToApply->eclipseCaseData(); + + RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS; + RimReservoirCellResultsStorage* gridCellResults = caseToApply->results(porosityModel); + RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); + + + std::vector fracDataVec; + + for (size_t fracCell : fracCells) + { + bool cellIsActive = activeCellInfo->isActive(fracCell); + + cvf::Vec3d localX; + cvf::Vec3d localY; + cvf::Vec3d localZ; + std::vector > planeCellPolygons; + bool isPlanIntersected = planeCellIntersectionPolygons(fracCell, planeCellPolygons, localX, localY, localZ); + if (!isPlanIntersected || planeCellPolygons.size() == 0) continue; + + //Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon/stimPlan mesh already is located) + cvf::Mat4f invertedTransMatrix = transformMatrix().getInverted(); + for (std::vector & planeCellPolygon : planeCellPolygons) + { + for (cvf::Vec3d& v : planeCellPolygon) + { + v.transformPoint(static_cast(invertedTransMatrix)); + } + } + + cvf::Vec3d localZinFracPlane; + localZinFracPlane = localZ; + localZinFracPlane.transformVector(static_cast(invertedTransMatrix)); + cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO; + directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1)); + directionOfLength.normalize(); + + RigFractureData fracData; + fracData.reservoirCellIndex = fracCell; + + + + std::vector fracPolygon = attachedFractureDefinition()->fracturePolygon(unitForExport); + std::vector > polygonsDescribingFractureInCell; + + double area; + std::vector areaOfFractureParts; + std::vector areaXvalueOfFractureParts; + + + for (std::vector planeCellPolygon : planeCellPolygons) + { + + for (int i = 0; i < stimPlanParameterValues.size(); i++) + { + double stimPlanParameterValue = stimPlanParameterValues[i]; + if (stimPlanParameterValue != 0) + { + std::vector stimPlanCell = stimPlanCellsAsPolygons[i]; + std::vector >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell, planeCellPolygon); + if (clippedStimPlanPolygons.size() > 0) + { + for (auto clippedStimPlanPolygon : clippedStimPlanPolygons) + { + area = cvf::GeometryTools::polygonAreaNormal3D(clippedStimPlanPolygon).length(); + areaOfFractureParts.push_back(area); + areaXvalueOfFractureParts.push_back(area*stimPlanParameterValue); + } + } + } + } + + + } + + if (areaXvalueOfFractureParts.size() > 0) + { + double fractureCellArea = 0.0; + for (double area : areaOfFractureParts) fractureCellArea += area; + double fractureCellAreaXvalue = 0.0; + for (double areaXvalue : areaXvalueOfFractureParts) fractureCellAreaXvalue += area; + double upscaledValue = fractureCellAreaXvalue / fractureCellArea; + + fracData.upscaledStimPlanValue = upscaledValue; + fracDataVec.push_back(fracData); + } + + + } + + m_rigFracture->setFractureData(fracDataVec); + + + +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ProjectDataModel/RimFracture.h b/ApplicationCode/ProjectDataModel/RimFracture.h index 51b4f1ef90..19e47a50fc 100644 --- a/ApplicationCode/ProjectDataModel/RimFracture.h +++ b/ApplicationCode/ProjectDataModel/RimFracture.h @@ -85,6 +85,8 @@ public: std::vector getPotentiallyFracturedCells(); void computeTransmissibility(RimEclipseCase* caseToApply); + void computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply); + virtual void fieldChangedByUi(const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue) override; cvf::Vec3d fracturePosition() const; diff --git a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp index 93d0c44e76..46821900a4 100644 --- a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp +++ b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp @@ -39,6 +39,8 @@ #include #include +#include +#include "RivWellFracturePartMgr.h" @@ -311,7 +313,6 @@ void RimStimPlanFractureTemplate::loadDataAndUpdate() updateConnectedEditors(); } - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -320,6 +321,24 @@ std::vector> RimStimPlanFractureTemplate::getDataAtTimeIndex return m_stimPlanFractureDefinitionData->getDataAtTimeIndex(resultName, unitName, timeStepIndex); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector> RimStimPlanFractureTemplate::getMirroredDataAtTimeIndex(const QString& resultName, const QString& unitName, size_t timeStepIndex) const +{ + std::vector> notMirrordedData = m_stimPlanFractureDefinitionData->getDataAtTimeIndex(resultName, unitName, timeStepIndex); + std::vector> mirroredData; + + for (std::vector depthData : notMirrordedData) + { + std::vector mirrordDepthData = RivWellFracturePartMgr::mirrorDataAtSingleDepth(depthData); + mirroredData.push_back(mirrordDepthData); + } + + + return mirroredData; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -709,6 +728,46 @@ void RimStimPlanFractureTemplate::computeMinMax(const QString& resultName, const } } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimStimPlanFractureTemplate::getStimPlanDataAsPolygonsAndValues(std::vector > &cellsAsPolygons, std::vector ¶meterValues, const QString& resultName, const QString& unitName, size_t timeStepIndex) +{ + + std::vector> propertyValuesAtTimeStep = getMirroredDataAtTimeIndex(resultName, unitName, timeStepIndex); + + cellsAsPolygons.clear(); + parameterValues.clear(); + + //TODO: Code partly copied from RivWellFracturePartMgr - can this be combined in some function? + std::vector depthCoordsAtNodes = adjustedDepthCoordsAroundWellPathPosition(); + std::vector xCoordsAtNodes = getNegAndPosXcoords(); + + //Cells are around nodes instead of between nodes + std::vector xCoords; + for (int i = 0; i < xCoordsAtNodes.size() - 1; i++) xCoords.push_back((xCoordsAtNodes[i] + xCoordsAtNodes[i + 1]) / 1); + std::vector depthCoords; + for (int i = 0; i < depthCoordsAtNodes.size() - 1; i++) depthCoords.push_back((depthCoordsAtNodes[i] + depthCoordsAtNodes[i + 1]) / 1); + + for (int i = 0; i < xCoords.size() - 1; i++) + { + for (int j = 0; j < depthCoords.size() - 1; j++) + { + std::vector cellAsPolygon; + cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i]), static_cast(depthCoords[j]), 0.0)); + cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i + 1]), static_cast(depthCoords[j]), 0.0)); + cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i + 1]), static_cast(depthCoords[j + 1]), 0.0)); + cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i]), static_cast(depthCoords[j + 1]), 0.0)); + cellsAsPolygons.push_back(cellAsPolygon); + //TODO: Values for both neg and pos x values... + parameterValues.push_back(propertyValuesAtTimeStep[j+1][i+1]); //TODO test that this value exsist... + + } + } + + +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h index 54448ea0b3..afb41286a4 100644 --- a/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h +++ b/ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h @@ -62,7 +62,6 @@ public: void fractureGeometry(std::vector* nodeCoords, std::vector* triangleIndices, caf::AppEnum< RimDefines::UnitSystem > fractureUnit) override; std::vector fracturePolygon(caf::AppEnum< RimDefines::UnitSystem > fractureUnit); - void sortPolygon(std::vector &polygon); std::vector getNegAndPosXcoords(); @@ -70,13 +69,13 @@ public: std::vector getStimPlanTimeValues(); std::vector > getStimPlanPropertyNamesUnits() const; void computeMinMax(const QString& resultName, const QString& unitName, double* minValue, double* maxValue) const; + void getStimPlanDataAsPolygonsAndValues(std::vector > &cellsAsPolygons, std::vector ¶meterValue, const QString& resultName, const QString& unitName, size_t timeStepIndex); void loadDataAndUpdate(); void setDefaultsBasedOnXMLfile(); - std::vector> getDataAtTimeIndex(const QString& resultName, const QString& unitName, size_t timeStepIndex) const; - - + std::vector> getMirroredDataAtTimeIndex(const QString& resultName, const QString& unitName, size_t timeStepIndex) const; + virtual QList calculateValueOptions(const caf::PdmFieldHandle* fieldNeedingOptions, bool* useOptionsOnly) override; protected: @@ -88,21 +87,19 @@ private: void updateUiTreeName(); void readStimPlanXMLFile(QString * errorMessage); - - size_t readStimplanGridAndTimesteps(QXmlStreamReader &xmlStream); - static double getAttributeValueDouble(QXmlStreamReader &xmlStream, QString parameterName); static QString getAttributeValueString(QXmlStreamReader &xmlStream, QString parameterName); void getGriddingValues(QXmlStreamReader &xmlStream, std::vector& gridValues, size_t& startNegValues); std::vector> getAllDepthDataAtTimeStep(QXmlStreamReader &xmlStream, size_t startingNegValuesXs); - caf::PdmField m_stimPlanFileName; - cvf::ref m_stimPlanFractureDefinitionData; - bool numberOfParameterValuesOK(std::vector> propertyValuesAtTimestep); bool setPropertyForPolygonDefault(); void setDepthOfWellPathAtFracture(); QString getUnitForStimPlanParameter(QString parameterName); + caf::PdmField m_stimPlanFileName; + cvf::ref m_stimPlanFractureDefinitionData; + + }; diff --git a/ApplicationCode/ReservoirDataModel/RigFracture.h b/ApplicationCode/ReservoirDataModel/RigFracture.h index 18ea4a5fbf..b1b8106d00 100644 --- a/ApplicationCode/ReservoirDataModel/RigFracture.h +++ b/ApplicationCode/ReservoirDataModel/RigFracture.h @@ -46,6 +46,7 @@ public: bool cellIsActive; + double upscaledStimPlanValue; };