///////////////////////////////////////////////////////////////////////////////// // // 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 "RigEclipseResultAddress.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 RimFracture* fracture ) : m_case( caseToApply ) , m_fractureTransform( fractureTransform ) , m_fractureSkinFactor( skinFactor ) , m_cDarcy( cDarcy ) , m_stimPlanCell( stimPlanCell ) , m_fracture( fracture ) { } void RigEclipseToStimPlanCellTransmissibilityCalculator::computeValues( const std::set& reservoirCellIndicesOpenForFlow ) { 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; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- 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 || std::isinf( m_stimPlanCell.getConductivityValue() ) ) return; const RigEclipseCaseData* eclipseCaseData = m_case->eclipseCaseData(); RiaDefines::PorosityModelType porosityModel = RiaDefines::PorosityModelType::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; std::vector stimPlanPolygon = m_stimPlanCell.getPolygon(); for ( const std::vector& planeCellPolygon : planeCellPolygons ) { std::vector> clippedPolygons = RigCellGeometryTools::intersectionWithPolygon( planeCellPolygon, stimPlanPolygon ); for ( const std::vector& clippedPolygon : clippedPolygons ) { polygonsForStimPlanCellInEclipseCell.push_back( clippedPolygon ); } } if ( polygonsForStimPlanCellInEclipseCell.empty() ) continue; std::vector areaOfFractureParts; std::vector lengthXareaOfFractureParts; double Ax = 0.0; double Ay = 0.0; double Az = 0.0; for ( const std::vector& fracturePartPolygon : polygonsForStimPlanCellInEclipseCell ) { cvf::Vec3d areaVector = cvf::GeometryTools::polygonAreaNormal3D( fracturePartPolygon ); double area = areaVector.length(); areaOfFractureParts.push_back( area ); double 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; // Guard against numerical issues for very small intersections if ( std::isnan( totalAreaXLength ) || std::isnan( fractureArea ) ) continue; 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->cell( 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 ); cvf::Vec3d transmissibilityVector( transmissibility_X, transmissibility_Y, transmissibility_Z ); transmissibility = calculateTransmissibility( transmissibilityVector, fractureArea ); 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 { const RigMainGrid* mainGrid = m_case->eclipseCaseData()->mainGrid(); if ( !mainGrid ) return {}; cvf::BoundingBox polygonBBox; for ( const cvf::Vec3d& nodeCoord : polygon ) { polygonBBox.add( nodeCoord ); } std::vector cellIndices = mainGrid->findIntersectingCells( polygonBBox ); std::vector cellIndicesToLeafCells; for ( const size_t& index : cellIndices ) { const RigCell& cell = mainGrid->cell( index ); if ( cell.isInvalid() ) continue; if ( !cell.subGrid() ) { cellIndicesToLeafCells.push_back( index ); } } return cellIndicesToLeafCells; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::ref RigEclipseToStimPlanCellTransmissibilityCalculator::createResultAccessor( const RimEclipseCase* eclipseCase, const QString& uiResultName ) { RiaDefines::PorosityModelType porosityModel = RiaDefines::PorosityModelType::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 ) ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- double RigEclipseToStimPlanCellTransmissibilityCalculator::calculateTransmissibility( const cvf::Vec3d& transmissibilityVector, double fractureArea ) { return transmissibilityVector.length(); }