///////////////////////////////////////////////////////////////////////////////// // // 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 "RigEclipseToStimPlanCellTransmissibilityCalculator.h" #include "RigActiveCellInfo.h" #include "RigCaseCellResultsData.h" #include "RigCellGeometryTools.h" #include "RigEclipseCaseData.h" #include "RigFractureCell.h" #include "RigFractureTransmissibilityEquations.h" #include "RigHexIntersectionTools.h" #include "RigMainGrid.h" #include "RigResultAccessorFactory.h" #include "RimEclipseCase.h" #include "RimFracture.h" #include "RiaLogging.h" #include "cvfGeometryTools.h" //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RigEclipseToStimPlanCellTransmissibilityCalculator::RigEclipseToStimPlanCellTransmissibilityCalculator( const RimEclipseCase* caseToApply, cvf::Mat4d fractureTransform, double skinFactor, double cDarcy, const RigFractureCell& stimPlanCell, const std::set& reservoirCellIndicesOpenForFlow, const RimFracture* fracture) : m_case(caseToApply) , m_fractureTransform(fractureTransform) , m_fractureSkinFactor(skinFactor) , m_cDarcy(cDarcy) , m_stimPlanCell(stimPlanCell) , m_fracture(fracture) { calculateStimPlanCellsMatrixTransmissibility(reservoirCellIndicesOpenForFlow); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& RigEclipseToStimPlanCellTransmissibilityCalculator::globalIndiciesToContributingEclipseCells() const { return m_globalIndiciesToContributingEclipseCells; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& RigEclipseToStimPlanCellTransmissibilityCalculator::contributingEclipseCellTransmissibilities() const { return m_contributingEclipseCellTransmissibilities; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& RigEclipseToStimPlanCellTransmissibilityCalculator::contributingEclipseCellIntersectionAreas() const { return m_contributingEclipseCellAreas; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& RigEclipseToStimPlanCellTransmissibilityCalculator::contributingEclipseCellPermeabilities() const { return m_contributingEclipseCellPermeabilities; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- double RigEclipseToStimPlanCellTransmissibilityCalculator::areaOpenForFlow() const { double area = 0.0; for (const auto& areaForOneEclipseCell : m_contributingEclipseCellAreas) { area += areaForOneEclipseCell; } return area; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- double RigEclipseToStimPlanCellTransmissibilityCalculator::aggregatedMatrixTransmissibility() const { double totalTransmissibility = 0.0; for (const auto& trans : m_contributingEclipseCellTransmissibilities) { totalTransmissibility += trans; } return totalTransmissibility; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const RigFractureCell& RigEclipseToStimPlanCellTransmissibilityCalculator::fractureCell() const { return m_stimPlanCell; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RigEclipseToStimPlanCellTransmissibilityCalculator::requiredResultNames() { std::vector resultNames; resultNames.push_back("PERMX"); resultNames.push_back("PERMY"); resultNames.push_back("PERMZ"); resultNames.push_back("DX"); resultNames.push_back("DY"); resultNames.push_back("DZ"); return resultNames; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RigEclipseToStimPlanCellTransmissibilityCalculator::optionalResultNames() { std::vector resultNames; resultNames.push_back("NTG"); return resultNames; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsMatrixTransmissibility( const std::set& reservoirCellIndicesOpenForFlow) { // Not calculating flow into fracture if stimPlan cell cond value is 0 (assumed to be outside the fracture): if (m_stimPlanCell.getConductivityValue() < 1e-7) return; const RigEclipseCaseData* eclipseCaseData = m_case->eclipseCaseData(); RiaDefines::PorosityModelType porosityModel = RiaDefines::MATRIX_MODEL; cvf::ref dataAccessObjectDx = createResultAccessor(m_case, "DX"); cvf::ref dataAccessObjectDy = createResultAccessor(m_case, "DY"); cvf::ref dataAccessObjectDz = createResultAccessor(m_case, "DZ"); if (dataAccessObjectDx.isNull() || dataAccessObjectDy.isNull() || dataAccessObjectDz.isNull()) { RiaLogging::error("Data for DX/DY/DZ is not complete, and these values are required for export of COMPDAT. Make sure " "'Preferences->Compute DEPTH Related Properties' is checked."); return; } cvf::ref dataAccessObjectPermX = createResultAccessor(m_case, "PERMX"); cvf::ref dataAccessObjectPermY = createResultAccessor(m_case, "PERMY"); cvf::ref dataAccessObjectPermZ = createResultAccessor(m_case, "PERMZ"); if (dataAccessObjectPermX.isNull() || dataAccessObjectPermY.isNull() || dataAccessObjectPermZ.isNull()) { RiaLogging::error("Data for PERMX/PERMY/PERMZ is not complete, and these values are required for export of COMPDAT."); return; } cvf::ref dataAccessObjectNTG = createResultAccessor(m_case, "NTG"); const RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); std::vector stimPlanPolygonTransformed; for (cvf::Vec3d v : m_stimPlanCell.getPolygon()) { v.transformPoint(m_fractureTransform); stimPlanPolygonTransformed.push_back(v); } std::vector reservoirCellIndices = getPotentiallyFracturedCellsForPolygon(stimPlanPolygonTransformed); for (size_t reservoirCellIndex : reservoirCellIndices) { const RigMainGrid* mainGrid = m_case->eclipseCaseData()->mainGrid(); if (!m_fracture->isEclipseCellOpenForFlow(mainGrid, reservoirCellIndicesOpenForFlow, reservoirCellIndex)) continue; std::array hexCorners; mainGrid->cellCornerVertices(reservoirCellIndex, hexCorners.data()); std::vector> planeCellPolygons; bool isPlanIntersected = RigHexIntersectionTools::planeHexIntersectionPolygons(hexCorners, m_fractureTransform, planeCellPolygons); if (!isPlanIntersected || planeCellPolygons.empty()) continue; cvf::Vec3d localX; cvf::Vec3d localY; cvf::Vec3d localZ; RigCellGeometryTools::findCellLocalXYZ(hexCorners, localX, localY, localZ); // Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon already is located) cvf::Mat4d invertedTransMatrix = m_fractureTransform.getInverted(); for (std::vector& planeCellPolygon : planeCellPolygons) { for (cvf::Vec3d& v : planeCellPolygon) { v.transformPoint(invertedTransMatrix); } } std::vector> polygonsForStimPlanCellInEclipseCell; cvf::Vec3d areaVector; std::vector stimPlanPolygon = m_stimPlanCell.getPolygon(); for (const std::vector& planeCellPolygon : planeCellPolygons) { std::vector> clippedPolygons = RigCellGeometryTools::intersectPolygons(planeCellPolygon, stimPlanPolygon); for (const std::vector& clippedPolygon : clippedPolygons) { polygonsForStimPlanCellInEclipseCell.push_back(clippedPolygon); } } if (polygonsForStimPlanCellInEclipseCell.empty()) continue; std::vector areaOfFractureParts; double length; std::vector lengthXareaOfFractureParts; double Ax = 0.0; double Ay = 0.0; double Az = 0.0; for (const std::vector& fracturePartPolygon : polygonsForStimPlanCellInEclipseCell) { areaVector = cvf::GeometryTools::polygonAreaNormal3D(fracturePartPolygon); double area = areaVector.length(); areaOfFractureParts.push_back(area); length = RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(fracturePartPolygon); lengthXareaOfFractureParts.push_back(length * area); cvf::Plane fracturePlane; fracturePlane.setFromPointAndNormal(static_cast(m_fractureTransform.translation()), static_cast(m_fractureTransform.col(2))); Ax += fabs(area * (fracturePlane.normal().dot(localY))); Ay += fabs(area * (fracturePlane.normal().dot(localX))); Az += fabs(area * (fracturePlane.normal().dot(localZ))); } double fractureArea = 0.0; for (double area : areaOfFractureParts) fractureArea += area; double totalAreaXLength = 0.0; for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea; double fractureAreaWeightedlength = totalAreaXLength / fractureArea; // Transmissibility for inactive cells is set to zero // Inactive cells must be include in order to compute the fractured area correctly double transmissibility = 0.0; bool isActive = true; { // Use main grid cell to evaluate if a cell is active or not. // All cells in temporary grids are active const RigCell& cell = mainGrid->globalCellArray()[reservoirCellIndex]; size_t mainGridReservoirIndex = cell.mainGridCellIndex(); if (!activeCellInfo->isActive(mainGridReservoirIndex)) { isActive = false; } } double matrixPermeability = 0.0; if (isActive) { double permX = dataAccessObjectPermX->cellScalarGlobIdx(reservoirCellIndex); double permY = dataAccessObjectPermY->cellScalarGlobIdx(reservoirCellIndex); double permZ = dataAccessObjectPermZ->cellScalarGlobIdx(reservoirCellIndex); double dx = dataAccessObjectDx->cellScalarGlobIdx(reservoirCellIndex); double dy = dataAccessObjectDy->cellScalarGlobIdx(reservoirCellIndex); double dz = dataAccessObjectDz->cellScalarGlobIdx(reservoirCellIndex); double NTG = 1.0; if (dataAccessObjectNTG.notNull()) { NTG = dataAccessObjectNTG->cellScalarGlobIdx(reservoirCellIndex); } double transmissibility_X = RigFractureTransmissibilityEquations::matrixToFractureTrans( permY, NTG, Ay, dx, m_fractureSkinFactor, fractureAreaWeightedlength, m_cDarcy); double transmissibility_Y = RigFractureTransmissibilityEquations::matrixToFractureTrans( permX, NTG, Ax, dy, m_fractureSkinFactor, fractureAreaWeightedlength, m_cDarcy); double transmissibility_Z = RigFractureTransmissibilityEquations::matrixToFractureTrans( permZ, 1.0, Az, dz, m_fractureSkinFactor, fractureAreaWeightedlength, m_cDarcy); transmissibility = sqrt(transmissibility_X * transmissibility_X + transmissibility_Y * transmissibility_Y + transmissibility_Z * transmissibility_Z); matrixPermeability = RigFractureTransmissibilityEquations::matrixPermeability(permX, permY, NTG); } m_globalIndiciesToContributingEclipseCells.push_back(reservoirCellIndex); m_contributingEclipseCellTransmissibilities.push_back(transmissibility); m_contributingEclipseCellAreas.push_back(fractureArea); m_contributingEclipseCellPermeabilities.push_back(matrixPermeability); } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RigEclipseToStimPlanCellTransmissibilityCalculator::getPotentiallyFracturedCellsForPolygon( const std::vector& polygon) const { std::vector cellIndices; const RigMainGrid* mainGrid = m_case->eclipseCaseData()->mainGrid(); if (!mainGrid) return cellIndices; cvf::BoundingBox polygonBBox; for (const cvf::Vec3d& nodeCoord : polygon) { polygonBBox.add(nodeCoord); } mainGrid->findIntersectingCells(polygonBBox, &cellIndices); std::vector cellIndicesToLeafCells; for (const size_t& index : cellIndices) { const RigCell& cell = mainGrid->globalCellArray()[index]; if (!cell.subGrid()) { cellIndicesToLeafCells.push_back(index); } } return cellIndicesToLeafCells; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::ref RigEclipseToStimPlanCellTransmissibilityCalculator::createResultAccessor(const RimEclipseCase* eclipseCase, const QString& uiResultName) { RiaDefines::PorosityModelType porosityModel = RiaDefines::MATRIX_MODEL; const RigEclipseCaseData* eclipseCaseData = eclipseCase->eclipseCaseData(); // Create result accessor object for main grid at time step zero (static result date is always at first time step return RigResultAccessorFactory::createFromResultAddress(eclipseCaseData, 0, porosityModel, 0, RigEclipseResultAddress(uiResultName)); }