From 56e84a43eb43566dcf269a56dca198a97064c67f Mon Sep 17 00:00:00 2001 From: astridkbjorke Date: Tue, 21 Mar 2017 13:43:44 +0100 Subject: [PATCH] pre-proto - Refactoring: Moiving functions related to transmissibility calculations from RimFracture to new class RigFractureTransCalc --- .../FileInterface/RifEclipseExportTools.cpp | 7 +- .../ProjectDataModel/RimFracture.cpp | 567 +---------------- .../ProjectDataModel/RimFracture.h | 14 +- .../ReservoirDataModel/CMakeLists_files.cmake | 2 + .../RigFractureTransCalc.cpp | 592 ++++++++++++++++++ .../ReservoirDataModel/RigFractureTransCalc.h | 86 +++ 6 files changed, 696 insertions(+), 572 deletions(-) create mode 100644 ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp create mode 100644 ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h diff --git a/ApplicationCode/FileInterface/RifEclipseExportTools.cpp b/ApplicationCode/FileInterface/RifEclipseExportTools.cpp index 26be1e16e8..547a2a47bf 100644 --- a/ApplicationCode/FileInterface/RifEclipseExportTools.cpp +++ b/ApplicationCode/FileInterface/RifEclipseExportTools.cpp @@ -39,6 +39,7 @@ #include #include #include +#include "RigFractureTransCalc.h" @@ -109,7 +110,8 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c for (RimFracture* fracture : fractures) { - fracture->computeTransmissibility(caseToApply); + //TODO: Check that there is a fracture template available for given fracture.... + RigFractureTransCalc::computeTransmissibility(caseToApply, fracture); std::vector fracDataVector = fracture->attachedRigFracture()->fractureData(); for (RigFractureData fracData : fracDataVector) @@ -124,7 +126,6 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c for (RimFracture* fracture : fractures) { RiaLogging::debug(QString("Writing COMPDAT values for fracture %1").arg(fracture->name())); - fracture->computeTransmissibility(caseToApply); std::vector fracDataVector = fracture->attachedRigFracture()->fractureData(); for (RigFractureData fracData : fracDataVector) @@ -163,7 +164,7 @@ void RifEclipseExportTools::performStimPlanUpscalingAndPrintResults(const std::v for (RimFracture* fracture : fractures) //For testing upscaling... { - fracture->computeUpscaledPropertyFromStimPlan(caseToApply, resultName, resultUnit, timeStepIndex); + RigFractureTransCalc::computeUpscaledPropertyFromStimPlan(caseToApply, fracture, resultName, resultUnit, timeStepIndex); std::vector fracDataVector = fracture->attachedRigFracture()->fractureData(); out << qSetFieldWidth(4); diff --git a/ApplicationCode/ProjectDataModel/RimFracture.cpp b/ApplicationCode/ProjectDataModel/RimFracture.cpp index 3eb257da2a..1855d8d560 100644 --- a/ApplicationCode/ProjectDataModel/RimFracture.cpp +++ b/ApplicationCode/ProjectDataModel/RimFracture.cpp @@ -132,6 +132,14 @@ const std::vector& RimFracture::nodeCoords() const return m_rigFracture->nodeCoords(); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimFracture::setFractureData(std::vector fractureDataVector) +{ + m_rigFracture->setFractureData(fractureDataVector); +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -276,565 +284,6 @@ cvf::Mat4f RimFracture::transformMatrix() return m; } - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -//TODO: Make static and move to another class -void RimFracture::computeTransmissibility(RimEclipseCase* caseToApply) -{ - //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: Handle case with finite conductivity in fracture - if (attachedFractureDefinition()) - { - if (attachedFractureDefinition()->fractureConductivity == RimFractureTemplate::FINITE_CONDUCTIVITY) - { - qDebug() << "Transimssibility for finite conductity in fracture not yet implemented."; - qDebug() << "Performing calculation for infinite conductivity instead."; - } - } - - - RigEclipseCaseData* eclipseCaseData = caseToApply->eclipseCaseData(); - - RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS; - RimReservoirCellResultsStorage* gridCellResults = caseToApply->results(porosityModel); - - size_t scalarSetIndex; - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DX"); - cvf::ref dataAccessObjectDx = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DX"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DY"); - cvf::ref dataAccessObjectDy = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DY"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DZ"); - cvf::ref dataAccessObjectDz = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DZ"); //assuming 0 time step and main grid (so grid index =0) - - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMX"); - cvf::ref dataAccessObjectPermX = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMX"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMY"); - cvf::ref dataAccessObjectPermY = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMY"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMZ"); - cvf::ref dataAccessObjectPermZ = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMZ"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "NTG"); - cvf::ref dataAccessObjectNTG = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "NTG"); //assuming 0 time step and main grid (so grid index =0) - - - RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); - - - std::vector fracDataVec; - - std::vector fracCells = getPotentiallyFracturedCells(); - - for (size_t fracCell : fracCells) - { - bool cellIsActive = activeCellInfo->isActive(fracCell); - - double permX = dataAccessObjectPermX->cellScalarGlobIdx(fracCell); - double permY = dataAccessObjectPermY->cellScalarGlobIdx(fracCell); - double permZ = dataAccessObjectPermZ->cellScalarGlobIdx(fracCell); - - double dx = dataAccessObjectDx->cellScalarGlobIdx(fracCell); - double dy = dataAccessObjectDy->cellScalarGlobIdx(fracCell); - double dz = dataAccessObjectDz->cellScalarGlobIdx(fracCell); - - double NTG = dataAccessObjectNTG->cellScalarGlobIdx(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 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; - - double transmissibility; - double fractureArea = 0.0; - double fractureAreaWeightedlength = 0.0; - double Ax = 0.0; - double Ay = 0.0; - double Az = 0.0; - double skinfactor = 0.0; - double transmissibility_X = 0.0; - double transmissibility_Y = 0.0; - double transmissibility_Z = 0.0; - - if (attachedFractureDefinition()) - { - std::vector fracPolygon = attachedFractureDefinition()->fracturePolygon(unitForExport); - - std::vector fracPolygonDouble; - for (auto v : fracPolygon) fracPolygonDouble.push_back(static_cast(v)); - - std::vector > polygonsDescribingFractureInCell; - cvf::Vec3d areaVector; - - for (std::vector planeCellPolygon : planeCellPolygons) - { - std::vector >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, fracPolygonDouble); - for (std::vector clippedPolygon : clippedPolygons) - { - polygonsDescribingFractureInCell.push_back(clippedPolygon); - } - } - - double area; - std::vector areaOfFractureParts; - double length; - std::vector lengthXareaOfFractureParts; - - for (std::vector fracturePartPolygon : polygonsDescribingFractureInCell) - { - areaVector = cvf::GeometryTools::polygonAreaNormal3D(fracturePartPolygon); - area = areaVector.length(); - areaOfFractureParts.push_back(area); - length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon); - lengthXareaOfFractureParts.push_back(length * area); - - cvf::Plane fracturePlane; - cvf::Mat4f m = transformMatrix(); - bool isCellIntersected = false; - - fracturePlane.setFromPointAndNormal(static_cast(m.translation()), - static_cast(m.col(2))); - - Ax += abs(area*(fracturePlane.normal().dot(localY))); - Ay += abs(area*(fracturePlane.normal().dot(localX))); - Az += abs(area*(fracturePlane.normal().dot(localZ))); - //TODO: resulting values have only been checked for vertical fracture... - } - - fractureArea = 0.0; - for (double area : areaOfFractureParts) fractureArea += area; - - double totalAreaXLength = 0.0; - for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea; - fractureAreaWeightedlength = totalAreaXLength / fractureArea; - - double c = cvf::UNDEFINED_DOUBLE; - if (unitForExport == RimDefines::UNITS_METRIC) c = 0.00852702; - if (unitForExport == RimDefines::UNITS_FIELD) c = 0.00112712; - // TODO: Use value from RimReservoirCellResultsStorage? - - skinfactor = attachedFractureDefinition()->skinFactor; - - double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D; - - transmissibility_X = 8 * c * ( permY * NTG ) * Ay / (dx + slDivPi); - transmissibility_Y = 8 * c * ( permX * NTG ) * Ax / (dy + slDivPi); - transmissibility_Z = 8 * c * permZ * Az / (dz + slDivPi); - - transmissibility = sqrt(transmissibility_X * transmissibility_X - + transmissibility_Y * transmissibility_Y - + transmissibility_Z * transmissibility_Z); - } - else - { - transmissibility = cvf::UNDEFINED_DOUBLE; - } - - - fracData.transmissibility = transmissibility; - fracData.transmissibilities = cvf::Vec3d(transmissibility_X, transmissibility_Y, transmissibility_Z); - - fracData.totalArea = fractureArea; - fracData.projectedAreas = cvf::Vec3d(Ax, Ay, Az); - fracData.fractureLenght = fractureAreaWeightedlength; - - fracData.cellSizes = cvf::Vec3d(dx, dy, dz); - fracData.permeabilities = cvf::Vec3d(permX, permY, permZ); - fracData.NTG = NTG; - fracData.skinFactor = skinfactor; - fracData.cellIsActive = cellIsActive; - - //Since we loop over all potentially fractured cells, we only keep FractureData for cells where fracture have an non-zero area. - if (fractureArea > 1e-5) - { - fracDataVec.push_back(fracData); - } - - } - - m_rigFracture->setFractureData(fracDataVec); -} - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RimFracture::computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply, QString resultName, QString resultUnit, size_t timeStepIndex) -{ - - //TODO: A lot of common code with function for calculating transmissibility... - - 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; - } - - 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) - { - - RigFractureData fracData; - fracData.reservoirCellIndex = fracCell; - - double upscaledAritmStimPlanValue = cvf::UNDEFINED_DOUBLE; - double upscaledHarmStimPlanValue = cvf::UNDEFINED_DOUBLE; - caf::AppEnum< RimDefines::UnitSystem > unitSystem = RimDefines::UNITS_METRIC; - computeUpscaledPropertyFromStimPlanForEclipseCell(upscaledAritmStimPlanValue, upscaledHarmStimPlanValue, caseToApply, resultName, resultUnit, timeStepIndex, unitSystem, fracCell); - - if (upscaledAritmStimPlanValue != cvf::UNDEFINED_DOUBLE) - { - fracData.upscaledAritmStimPlanValue = upscaledAritmStimPlanValue; - fracData.upscaledHarmStimPlanValue = upscaledHarmStimPlanValue; - - fracDataVec.push_back(fracData); - } - } - - m_rigFracture->setFractureData(fracDataVec); - - - -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RimFracture::computeUpscaledPropertyFromStimPlanForEclipseCell(double &upscaledAritmStimPlanValue, double &upscaledHarmStimPlanValue, RimEclipseCase* caseToApply, QString resultName, QString resultUnit, size_t timeStepIndex, caf::AppEnum< RimDefines::UnitSystem > unitSystem, size_t cellIndex) -{ - //TODO: A lot of common code with function for calculating transmissibility... - - if (!attachedFractureDefinition()) return; - - RimStimPlanFractureTemplate* fracTemplateStimPlan; - if (dynamic_cast(attachedFractureDefinition())) - { - fracTemplateStimPlan = dynamic_cast(attachedFractureDefinition()); - } - else return; - - //TODO: UNITS! - - 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); - - - bool cellIsActive = activeCellInfo->isActive(cellIndex); - - cvf::Vec3d localX; - cvf::Vec3d localY; - cvf::Vec3d localZ; - std::vector > planeCellPolygons; - bool isPlanIntersected = planeCellIntersectionPolygons(cellIndex, planeCellPolygons, localX, localY, localZ); - if (!isPlanIntersected || planeCellPolygons.size() == 0) return; - - //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(); - - std::vector fracPolygon = attachedFractureDefinition()->fracturePolygon(unitSystem); - std::vector > polygonsDescribingFractureInCell; - - double area; - std::vector areaOfFractureParts; - std::vector valuesForFractureParts; - - 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); - valuesForFractureParts.push_back(stimPlanParameterValue); - } - } - } - } - } - - if (areaOfFractureParts.size() > 0) - { - upscaledAritmStimPlanValue = areaWeightedArithmeticAverage(areaOfFractureParts, valuesForFractureParts); - upscaledHarmStimPlanValue = areaWeightedHarmonicAverage(areaOfFractureParts, valuesForFractureParts); - } -} - - - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RimFracture::computeFlowInFracture(RimEclipseCase* caseToApply) -{ - - //TODO: A lot of common code with function for calculating transmissibility... - - if (!attachedFractureDefinition()) return; - - RimStimPlanFractureTemplate* fracTemplateStimPlan; - RimEllipseFractureTemplate* fracTemplateEllipse; - if (dynamic_cast(attachedFractureDefinition())) - { - fracTemplateStimPlan = dynamic_cast(attachedFractureDefinition()); - } - else if (dynamic_cast(attachedFractureDefinition())) - { - fracTemplateEllipse = 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: 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) - { - double K; //Conductivity - double w; //Fracture width - - if (fracTemplateEllipse) - { - //TODO: UNit handling... - K = fracTemplateEllipse->fractureConductivity(); - w = fracTemplateEllipse->width(); - } - else if (fracTemplateStimPlan) - { - QString resultName = "CONDUCTIVITY"; - QString resultUnit = "md-m"; - size_t timeStepIndex = 0; //TODO... - double upscaledAritmStimPlanValue = cvf::UNDEFINED_DOUBLE; - double upscaledHarmStimPlanValue = cvf::UNDEFINED_DOUBLE; - caf::AppEnum< RimDefines::UnitSystem > unitSystem = RimDefines::UNITS_METRIC; - computeUpscaledPropertyFromStimPlanForEclipseCell(upscaledAritmStimPlanValue, upscaledHarmStimPlanValue, caseToApply, resultName, resultUnit, timeStepIndex, unitSystem, fracCell); - K = (upscaledAritmStimPlanValue + upscaledHarmStimPlanValue) / 2; - - resultName = "WIDTH"; - resultUnit = "cm"; //TODO handle mm and cm! - computeUpscaledPropertyFromStimPlanForEclipseCell(upscaledAritmStimPlanValue, upscaledHarmStimPlanValue, caseToApply, resultName, resultUnit, timeStepIndex, unitSystem, fracCell); - w = (upscaledAritmStimPlanValue + upscaledHarmStimPlanValue) / 2; - - } - - } - - m_rigFracture->setFractureData(fracDataVec); - -} - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RimFracture::computeFlowIntoTransverseWell(RimEclipseCase* caseToApply) -{ - - //TODO: A lot of common code with function for calculating transmissibility... - - if (!attachedFractureDefinition()) return; - if (attachedFractureDefinition()->orientation == RimFractureTemplate::ALONG_WELL_PATH) return; - - double wellRadius = cvf::UNDEFINED_DOUBLE; - if (attachedFractureDefinition()->orientation == RimFractureTemplate::TRANSVERSE_WELL_PATH) - { - wellRadius = 0.0;//TODO read this value... - } - if (attachedFractureDefinition()->orientation == RimFractureTemplate::AZIMUTH) - { - - - } - - - -} - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimFracture::areaWeightedHarmonicAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts) -{ - //TODO: Unit test? - double fractureCellArea = 0.0; - for (double area : areaOfFractureParts) fractureCellArea += area; - - if (areaOfFractureParts.size() != valuesForFractureParts.size()) return cvf::UNDEFINED_DOUBLE; - - double fractureCellAreaDivvalue = 0.0; - for (int i = 0; i < valuesForFractureParts.size(); i++) - { - fractureCellAreaDivvalue += areaOfFractureParts[i] / valuesForFractureParts[i]; - } - - double upscaledValueHarmonic = fractureCellArea / fractureCellAreaDivvalue; - return upscaledValueHarmonic; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimFracture::areaWeightedArithmeticAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts) -{ - //TODO: Unit test? - double fractureCellArea = 0.0; - for (double area : areaOfFractureParts) fractureCellArea += area; - - if (areaOfFractureParts.size() != valuesForFractureParts.size()) return cvf::UNDEFINED_DOUBLE; - - double fractureCellAreaXvalue = 0.0; - for (int i = 0; i < valuesForFractureParts.size(); i++) - { - fractureCellAreaXvalue += areaOfFractureParts[i] * valuesForFractureParts[i]; - } - - double upscaledValueArithmetic = fractureCellAreaXvalue / fractureCellArea; - return upscaledValueArithmetic; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ProjectDataModel/RimFracture.h b/ApplicationCode/ProjectDataModel/RimFracture.h index 81ab9dcaed..29e2d50c8b 100644 --- a/ApplicationCode/ProjectDataModel/RimFracture.h +++ b/ApplicationCode/ProjectDataModel/RimFracture.h @@ -39,6 +39,7 @@ class RimEclipseCase; class RimEllipseFractureTemplate; class RivWellFracturePartMgr; class RimFractureTemplate; +class RigFractureData; //================================================================================================== /// @@ -81,18 +82,13 @@ public: const std::vector& triangleIndices() const; const std::vector& nodeCoords() const; + void setFractureData(std::vector fractureDataVector); std::vector getPotentiallyFracturedCells(); - void computeTransmissibility(RimEclipseCase* caseToApply); - void computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply, QString resultName, QString resultUnit, size_t timeStepIndex); - void computeUpscaledPropertyFromStimPlanForEclipseCell(double &upscaledAritmStimPlanValue, double &upscaledHarmStimPlanValue, RimEclipseCase* caseToApply, QString resultName, QString resultUnit, size_t timeStepIndex, caf::AppEnum< RimDefines::UnitSystem > unitSystem, size_t cellIndex); - void computeFlowInFracture(RimEclipseCase* caseToApply); - void computeFlowIntoTransverseWell(RimEclipseCase* caseToApply); - - double areaWeightedHarmonicAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts); - double areaWeightedArithmeticAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts); + //TODO: Move this function? + bool planeCellIntersectionPolygons(size_t cellindex, std::vector > & polygons, cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ); virtual void fieldChangedByUi(const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue) override; cvf::Vec3d fracturePosition() const; @@ -114,8 +110,6 @@ private: QString createOneBasedIJK() const; - //Functions for area calculations - should these be in separate class - bool planeCellIntersectionPolygons(size_t cellindex, std::vector > & polygons, cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ); protected: caf::PdmPtrField m_fractureTemplate; diff --git a/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake b/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake index 0a1bbc7cc3..b4b3782666 100644 --- a/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake +++ b/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake @@ -52,6 +52,7 @@ ${CEE_CURRENT_LIST_DIR}RigSummaryCaseData.h ${CEE_CURRENT_LIST_DIR}RigLasFileExporter.h ${CEE_CURRENT_LIST_DIR}RigSimulationWellCoordsAndMD.h ${CEE_CURRENT_LIST_DIR}RigFracture.h +${CEE_CURRENT_LIST_DIR}RigFractureTransCalc.h ${CEE_CURRENT_LIST_DIR}RigTesselatorTools.h ${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.h ${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.h @@ -103,6 +104,7 @@ ${CEE_CURRENT_LIST_DIR}RigSummaryCaseData.cpp ${CEE_CURRENT_LIST_DIR}RigLasFileExporter.cpp ${CEE_CURRENT_LIST_DIR}RigSimulationWellCoordsAndMD.cpp ${CEE_CURRENT_LIST_DIR}RigFracture.cpp +${CEE_CURRENT_LIST_DIR}RigFractureTransCalc.cpp ${CEE_CURRENT_LIST_DIR}RigTesselatorTools.cpp ${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.cpp ${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.cpp diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp new file mode 100644 index 0000000000..dea88c6323 --- /dev/null +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp @@ -0,0 +1,592 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2017 Statoil ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RigFractureTransCalc.h" +#include "RimFractureTemplate.h" +#include "RigEclipseCaseData.h" +#include "RimEclipseCase.h" +#include "RiaLogging.h" +#include "QString" +#include "RimReservoirCellResultsStorage.h" +#include "RigResultAccessorFactory.h" +#include "RimFracture.h" +#include "RigFracture.h" +#include "cvfGeometryTools.h" +#include "RigCellGeometryTools.h" +#include "RigActiveCellInfo.h" +#include "RimStimPlanFractureTemplate.h" + +#include +#include "RimEllipseFractureTemplate.h" +#include "cafAppEnum.h" + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +// RigFractureData::RigFractureData() +// { +// +// } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFractureTransCalc::RigFractureTransCalc() +{ +} + + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +//TODO: Make static and move to another class +void RigFractureTransCalc::computeTransmissibility(RimEclipseCase* caseToApply, RimFracture* fracture) +{ + //Get correct unit system: + RigEclipseCaseData::UnitsType caseUnit = caseToApply->eclipseCaseData()->unitsType(); + RimDefines::UnitSystem unitForExport; + + if (caseUnit == RigEclipseCaseData::UNITS_METRIC) + { + RiaLogging::debug(QString("Calculating transmissibilities in metric units")); + unitForExport = RimDefines::UNITS_METRIC; + } + else if (caseUnit == RigEclipseCaseData::UNITS_FIELD) + { + RiaLogging::debug(QString("Calculating transmissibilities in field units")); + unitForExport = RimDefines::UNITS_FIELD; + } + else + { + RiaLogging::error(QString("Unit system for case not supported for fracture export.")); + return; + } + + if (fracture->attachedFractureDefinition()->fractureConductivity == RimFractureTemplate::FINITE_CONDUCTIVITY) + { + RiaLogging::warning(QString("Transimssibility for finite conductity in fracture not yet implemented.")); + RiaLogging::warning(QString("Performing calculation for infinite conductivity instead.")); + } + + + RigEclipseCaseData* eclipseCaseData = caseToApply->eclipseCaseData(); + + RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS; + RimReservoirCellResultsStorage* gridCellResults = caseToApply->results(porosityModel); + + size_t scalarSetIndex; + scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DX"); + cvf::ref dataAccessObjectDx = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DX"); //assuming 0 time step and main grid (so grid index =0) + scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DY"); + cvf::ref dataAccessObjectDy = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DY"); //assuming 0 time step and main grid (so grid index =0) + scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DZ"); + cvf::ref dataAccessObjectDz = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DZ"); //assuming 0 time step and main grid (so grid index =0) + + scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMX"); + cvf::ref dataAccessObjectPermX = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMX"); //assuming 0 time step and main grid (so grid index =0) + scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMY"); + cvf::ref dataAccessObjectPermY = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMY"); //assuming 0 time step and main grid (so grid index =0) + scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMZ"); + cvf::ref dataAccessObjectPermZ = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMZ"); //assuming 0 time step and main grid (so grid index =0) + scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "NTG"); + cvf::ref dataAccessObjectNTG = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "NTG"); //assuming 0 time step and main grid (so grid index =0) + + RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); + + std::vector fracDataVec; + std::vector fracCells = fracture->getPotentiallyFracturedCells(); + + for (size_t fracCell : fracCells) + { + bool cellIsActive = activeCellInfo->isActive(fracCell); + + double permX = dataAccessObjectPermX->cellScalarGlobIdx(fracCell); + double permY = dataAccessObjectPermY->cellScalarGlobIdx(fracCell); + double permZ = dataAccessObjectPermZ->cellScalarGlobIdx(fracCell); + + double dx = dataAccessObjectDx->cellScalarGlobIdx(fracCell); + double dy = dataAccessObjectDy->cellScalarGlobIdx(fracCell); + double dz = dataAccessObjectDz->cellScalarGlobIdx(fracCell); + + double NTG = dataAccessObjectNTG->cellScalarGlobIdx(fracCell); + + cvf::Vec3d localX; + cvf::Vec3d localY; + cvf::Vec3d localZ; + std::vector > planeCellPolygons; + bool isPlanIntersected = fracture->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 already is located) + cvf::Mat4f invertedTransMatrix = fracture->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; + + double transmissibility; + double fractureArea = 0.0; + double fractureAreaWeightedlength = 0.0; + double Ax = 0.0; + double Ay = 0.0; + double Az = 0.0; + double skinfactor = 0.0; + double transmissibility_X = 0.0; + double transmissibility_Y = 0.0; + double transmissibility_Z = 0.0; + + + std::vector fracPolygon = fracture->attachedFractureDefinition()->fracturePolygon(unitForExport); + + std::vector fracPolygonDouble; + for (auto v : fracPolygon) fracPolygonDouble.push_back(static_cast(v)); + + std::vector > polygonsDescribingFractureInCell; + cvf::Vec3d areaVector; + + for (std::vector planeCellPolygon : planeCellPolygons) + { + std::vector >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, fracPolygonDouble); + for (std::vector clippedPolygon : clippedPolygons) + { + polygonsDescribingFractureInCell.push_back(clippedPolygon); + } + } + + double area; + std::vector areaOfFractureParts; + double length; + std::vector lengthXareaOfFractureParts; + + for (std::vector fracturePartPolygon : polygonsDescribingFractureInCell) + { + areaVector = cvf::GeometryTools::polygonAreaNormal3D(fracturePartPolygon); + area = areaVector.length(); + areaOfFractureParts.push_back(area); + length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon); + lengthXareaOfFractureParts.push_back(length * area); + + cvf::Plane fracturePlane; + cvf::Mat4f m = fracture->transformMatrix(); + bool isCellIntersected = false; + + fracturePlane.setFromPointAndNormal(static_cast(m.translation()), + static_cast(m.col(2))); + + Ax += abs(area*(fracturePlane.normal().dot(localY))); + Ay += abs(area*(fracturePlane.normal().dot(localX))); + Az += abs(area*(fracturePlane.normal().dot(localZ))); + //TODO: resulting values have only been checked for vertical fracture... + } + + fractureArea = 0.0; + for (double area : areaOfFractureParts) fractureArea += area; + + double totalAreaXLength = 0.0; + for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea; + fractureAreaWeightedlength = totalAreaXLength / fractureArea; + + double c = cvf::UNDEFINED_DOUBLE; + if (unitForExport == RimDefines::UNITS_METRIC) c = 0.00852702; + if (unitForExport == RimDefines::UNITS_FIELD) c = 0.00112712; + // TODO: Use value from RimReservoirCellResultsStorage? + + skinfactor = fracture->attachedFractureDefinition()->skinFactor; + + double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D; + + transmissibility_X = 8 * c * (permY * NTG) * Ay / (dx + slDivPi); + transmissibility_Y = 8 * c * (permX * NTG) * Ax / (dy + slDivPi); + transmissibility_Z = 8 * c * permZ * Az / (dz + slDivPi); + + transmissibility = sqrt(transmissibility_X * transmissibility_X + + transmissibility_Y * transmissibility_Y + + transmissibility_Z * transmissibility_Z); + + fracData.transmissibility = transmissibility; + fracData.transmissibilities = cvf::Vec3d(transmissibility_X, transmissibility_Y, transmissibility_Z); + + fracData.totalArea = fractureArea; + fracData.projectedAreas = cvf::Vec3d(Ax, Ay, Az); + fracData.fractureLenght = fractureAreaWeightedlength; + + fracData.cellSizes = cvf::Vec3d(dx, dy, dz); + fracData.permeabilities = cvf::Vec3d(permX, permY, permZ); + fracData.NTG = NTG; + fracData.skinFactor = skinfactor; + fracData.cellIsActive = cellIsActive; + + //Since we loop over all potentially fractured cells, we only keep FractureData for cells where fracture have an non-zero area. + if (fractureArea > 1e-5) + { + fracDataVec.push_back(fracData); + } + + } + + fracture->setFractureData(fracDataVec); +} + + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigFractureTransCalc::computeUpscaledPropertyFromStimPlanForEclipseCell(double &upscaledAritmStimPlanValue, double &upscaledHarmStimPlanValue, RimFracture* fracture, RimEclipseCase* caseToApply, QString resultName, QString resultUnit, size_t timeStepIndex, caf::AppEnum< RimDefines::UnitSystem > unitSystem, size_t cellIndex) +{ + //TODO: A lot of common code with function for calculating transmissibility... + + RimStimPlanFractureTemplate* fracTemplateStimPlan; + if (dynamic_cast(fracture->attachedFractureDefinition())) + { + fracTemplateStimPlan = dynamic_cast(fracture->attachedFractureDefinition()); + } + else return; + + //TODO: UNITS! + + 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 = fracture->getPotentiallyFracturedCells(); + + + RigEclipseCaseData* eclipseCaseData = caseToApply->eclipseCaseData(); + + RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS; + RimReservoirCellResultsStorage* gridCellResults = caseToApply->results(porosityModel); + RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); + + + bool cellIsActive = activeCellInfo->isActive(cellIndex); + + cvf::Vec3d localX; + cvf::Vec3d localY; + cvf::Vec3d localZ; + std::vector > planeCellPolygons; + bool isPlanIntersected = fracture->planeCellIntersectionPolygons(cellIndex, planeCellPolygons, localX, localY, localZ); + if (!isPlanIntersected || planeCellPolygons.size() == 0) return; + + //Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon/stimPlan mesh already is located) + cvf::Mat4f invertedTransMatrix = fracture->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(); + + std::vector fracPolygon = fracture->attachedFractureDefinition()->fracturePolygon(unitSystem); + std::vector > polygonsDescribingFractureInCell; + + double area; + std::vector areaOfFractureParts; + std::vector valuesForFractureParts; + + 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); + valuesForFractureParts.push_back(stimPlanParameterValue); + } + } + } + } + } + + if (areaOfFractureParts.size() > 0) + { + upscaledAritmStimPlanValue = areaWeightedArithmeticAverage(areaOfFractureParts, valuesForFractureParts); + upscaledHarmStimPlanValue = areaWeightedHarmonicAverage(areaOfFractureParts, valuesForFractureParts); + } +} + + + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigFractureTransCalc::areaWeightedHarmonicAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts) +{ + //TODO: Unit test? + double fractureCellArea = 0.0; + for (double area : areaOfFractureParts) fractureCellArea += area; + + if (areaOfFractureParts.size() != valuesForFractureParts.size()) return cvf::UNDEFINED_DOUBLE; + + double fractureCellAreaDivvalue = 0.0; + for (int i = 0; i < valuesForFractureParts.size(); i++) + { + fractureCellAreaDivvalue += areaOfFractureParts[i] / valuesForFractureParts[i]; + } + + double upscaledValueHarmonic = fractureCellArea / fractureCellAreaDivvalue; + return upscaledValueHarmonic; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigFractureTransCalc::areaWeightedArithmeticAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts) +{ + //TODO: Unit test? + double fractureCellArea = 0.0; + for (double area : areaOfFractureParts) fractureCellArea += area; + + if (areaOfFractureParts.size() != valuesForFractureParts.size()) return cvf::UNDEFINED_DOUBLE; + + double fractureCellAreaXvalue = 0.0; + for (int i = 0; i < valuesForFractureParts.size(); i++) + { + fractureCellAreaXvalue += areaOfFractureParts[i] * valuesForFractureParts[i]; + } + + double upscaledValueArithmetic = fractureCellAreaXvalue / fractureCellArea; + return upscaledValueArithmetic; +} + + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigFractureTransCalc::computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply, RimFracture* fracture, QString resultName, QString resultUnit, size_t timeStepIndex) +{ + + //TODO: A lot of common code with function for calculating transmissibility... + + RimStimPlanFractureTemplate* fracTemplateStimPlan; + if (dynamic_cast(fracture->attachedFractureDefinition())) + { + fracTemplateStimPlan = dynamic_cast(fracture->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 upscaled stimPlan values in metric units")); + unitForExport = RimDefines::UNITS_METRIC; + } + else if (caseUnit == RigEclipseCaseData::UNITS_FIELD) + { + RiaLogging::debug(QString("Calculating upscaled stimPlan values in field units")); + unitForExport = RimDefines::UNITS_FIELD; + } + else + { + RiaLogging::error(QString("Unit system for case not supported for fracture export.")); + return; + } + + 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 = fracture->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) + { + + RigFractureData fracData; + fracData.reservoirCellIndex = fracCell; + + double upscaledAritmStimPlanValue = cvf::UNDEFINED_DOUBLE; + double upscaledHarmStimPlanValue = cvf::UNDEFINED_DOUBLE; + caf::AppEnum< RimDefines::UnitSystem > unitSystem = RimDefines::UNITS_METRIC; + computeUpscaledPropertyFromStimPlanForEclipseCell(upscaledAritmStimPlanValue, upscaledHarmStimPlanValue, fracture, caseToApply, resultName, resultUnit, timeStepIndex, unitSystem, fracCell); + + if (upscaledAritmStimPlanValue != cvf::UNDEFINED_DOUBLE) + { + fracData.upscaledAritmStimPlanValue = upscaledAritmStimPlanValue; + fracData.upscaledHarmStimPlanValue = upscaledHarmStimPlanValue; + + fracDataVec.push_back(fracData); + } + } + + fracture->setFractureData(fracDataVec); + + + +} + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigFractureTransCalc::computeFlowInFracture(RimEclipseCase* caseToApply, RimFracture* fracture) +{ + + //TODO: A lot of common code with function for calculating transmissibility... + + RimStimPlanFractureTemplate* fracTemplateStimPlan; + RimEllipseFractureTemplate* fracTemplateEllipse; + if (dynamic_cast(fracture->attachedFractureDefinition())) + { + fracTemplateStimPlan = dynamic_cast(fracture->attachedFractureDefinition()); + } + else if (dynamic_cast(fracture->attachedFractureDefinition())) + { + fracTemplateEllipse = dynamic_cast(fracture->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 flow in fracture in metric units")); + unitForExport = RimDefines::UNITS_METRIC; + } + else if (caseUnit == RigEclipseCaseData::UNITS_FIELD) + { + RiaLogging::debug(QString("Calculating flow in fracture in field units")); + unitForExport = RimDefines::UNITS_FIELD; + } + else + { + RiaLogging::error(QString("Unit system for case not supported for fracture export.")); + return; + } + + + //TODO: A lot of common code with function above... Can be cleaned up...? + std::vector fracCells = fracture->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) + { + double K; //Conductivity + double w; //Fracture width + + if (fracTemplateEllipse) + { + //TODO: UNit handling... + K = fracTemplateEllipse->fractureConductivity(); + w = fracTemplateEllipse->width(); + } + else if (fracTemplateStimPlan) + { + QString resultName = "CONDUCTIVITY"; + QString resultUnit = "md-m"; + size_t timeStepIndex = 0; //TODO... + double upscaledAritmStimPlanValue = cvf::UNDEFINED_DOUBLE; + double upscaledHarmStimPlanValue = cvf::UNDEFINED_DOUBLE; + caf::AppEnum< RimDefines::UnitSystem > unitSystem = RimDefines::UNITS_METRIC; + computeUpscaledPropertyFromStimPlanForEclipseCell(upscaledAritmStimPlanValue, upscaledHarmStimPlanValue, fracture, caseToApply, resultName, resultUnit, timeStepIndex, unitSystem, fracCell); + K = (upscaledAritmStimPlanValue + upscaledHarmStimPlanValue) / 2; + + resultName = "WIDTH"; + resultUnit = "cm"; //TODO handle mm and cm! + computeUpscaledPropertyFromStimPlanForEclipseCell(upscaledAritmStimPlanValue, upscaledHarmStimPlanValue, fracture, caseToApply, resultName, resultUnit, timeStepIndex, unitSystem, fracCell); + w = (upscaledAritmStimPlanValue + upscaledHarmStimPlanValue) / 2; + + } + + } + + fracture->setFractureData(fracDataVec); + +} + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigFractureTransCalc::computeFlowIntoTransverseWell(RimEclipseCase* caseToApply, RimFracture* fracture) +{ + + //TODO: A lot of common code with function for calculating transmissibility... + + if (fracture->attachedFractureDefinition()->orientation == RimFractureTemplate::ALONG_WELL_PATH) return; + + double wellRadius = cvf::UNDEFINED_DOUBLE; + if (fracture->attachedFractureDefinition()->orientation == RimFractureTemplate::TRANSVERSE_WELL_PATH) + { + wellRadius = 0.0;//TODO read this value... + } + if (fracture->attachedFractureDefinition()->orientation == RimFractureTemplate::AZIMUTH) + { + + + } + + + +} + diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h new file mode 100644 index 0000000000..fb600ca25a --- /dev/null +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h @@ -0,0 +1,86 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2017 Statoil ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + + +#include "cvfBase.h" +#include "cvfObject.h" +#include "cvfMath.h" +#include "cvfVector3.h" + +#include +#include +#include "cafAppEnum.h" +#include "RimDefines.h" + + +class RimFracture; +class RimEclipseCase; + +// class RigFractureData +// { +// public: +// RigFractureData(); +// +// size_t reservoirCellIndex; +// double transmissibility; +// cvf::Vec3d transmissibilities; +// +// double totalArea; +// double fractureLenght; +// cvf::Vec3d projectedAreas; +// +// cvf::Vec3d permeabilities; +// cvf::Vec3d cellSizes; +// double NTG; +// double skinFactor; +// +// bool cellIsActive; +// +// //TODO: Used for upscaling - should be moved? +// double upscaledAritmStimPlanValue; +// double upscaledHarmStimPlanValue; +// +// +// }; + +//================================================================================================== +/// +//================================================================================================== +class RigFractureTransCalc : public cvf::Object +{ +public: + RigFractureTransCalc(); + + static void computeTransmissibility(RimEclipseCase* caseToApply, RimFracture* fracture); + + static void computeUpscaledPropertyFromStimPlanForEclipseCell(double &upscaledAritmStimPlanValue, double &upscaledHarmStimPlanValue, RimFracture* fracture, RimEclipseCase* caseToApply, QString resultName, QString resultUnit, size_t timeStepIndex, caf::AppEnum< RimDefines::UnitSystem > unitSystem, size_t cellIndex); + static double areaWeightedHarmonicAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts); + static double areaWeightedArithmeticAverage(std::vector areaOfFractureParts, std::vector valuesForFractureParts); + + static void computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply, RimFracture* fracture, QString resultName, QString resultUnit, size_t timeStepIndex); + + static void computeFlowInFracture(RimEclipseCase* caseToApply, RimFracture* fracture); + static void computeFlowIntoTransverseWell(RimEclipseCase* caseToApply, RimFracture* fracture); + + +private: + +}; +