///////////////////////////////////////////////////////////////////////////////// // // Copyright (C) 2023 Equinor 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 "RimFaultReactivationDataAccessorWellLogExtraction.h" #include "RiaEclipseUnitTools.h" #include "RiaInterpolationTools.h" #include "RiaLogging.h" #include "RigEclipseWellLogExtractor.h" #include "RigFaultReactivationModel.h" #include "RigFemAddressDefines.h" #include "RigFemPartCollection.h" #include "RigFemPartResultsCollection.h" #include "RigFemResultAddress.h" #include "RigFemScalarResultFrames.h" #include "RigGeoMechCaseData.h" #include "RigGeoMechWellLogExtractor.h" #include "RigGriddedPart3d.h" #include "RigResultAccessorFactory.h" #include "RigWellPath.h" #include "RimFaultReactivationEnums.h" #include "RimFracture.h" #include "RimGeoMechCase.h" #include "RimWellIADataAccess.h" #include "cvfVector3.h" #include #include //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RimFaultReactivationDataAccessorWellLogExtraction::RimFaultReactivationDataAccessorWellLogExtraction() { } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RimFaultReactivationDataAccessorWellLogExtraction::~RimFaultReactivationDataAccessorWellLogExtraction() { } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::pair RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( const std::vector& intersections, std::vector& values, const cvf::Vec3d& position, double gradient ) { // Fill in missing values fillInMissingValuesWithGradient( intersections, values, gradient ); return findValueAndPosition( intersections, values, position ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::pair RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( const std::vector& intersections, std::vector& values, const cvf::Vec3d& position, double seabedTemperature ) { // Fill in missing values fillInMissingValuesWithTopValue( intersections, values, seabedTemperature ); return findValueAndPosition( intersections, values, position ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::pair RimFaultReactivationDataAccessorWellLogExtraction::findValueAndPosition( const std::vector& intersections, const std::vector& values, const cvf::Vec3d& position ) { // Linear interpolation between two points auto lerp = []( const cvf::Vec3d& start, const cvf::Vec3d& end, double t ) { return start + t * ( end - start ); }; auto [topIdx, bottomIdx] = findIntersectionsForTvd( intersections, position.z() ); if ( topIdx != -1 && bottomIdx != -1 ) { double topValue = values[topIdx]; double bottomValue = values[bottomIdx]; if ( !std::isinf( topValue ) && !std::isinf( bottomValue ) ) { // Interpolate value from the two closest points. std::vector xs = { intersections[bottomIdx].z(), intersections[topIdx].z() }; std::vector ys = { values[bottomIdx], values[topIdx] }; double porBar = RiaInterpolationTools::linear( xs, ys, position.z() ); // Interpolate position from depth double fraction = RiaInterpolationTools::linear( xs, { 0.0, 1.0 }, position.z() ); cvf::Vec3d extractionPosition = lerp( intersections[bottomIdx], intersections[topIdx], fraction ); return { porBar, extractionPosition }; } } else if ( position.z() <= intersections.back().z() ) { return { values.back(), intersections.back() }; } return { std::numeric_limits::infinity(), cvf::Vec3d::UNDEFINED }; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::pair RimFaultReactivationDataAccessorWellLogExtraction::findIntersectionsForTvd( const std::vector& intersections, double tvd ) { int topIdx = -1; int bottomIdx = -1; if ( intersections.size() >= 2 ) { for ( size_t i = 1; i < intersections.size(); i++ ) { auto top = intersections[i - 1]; auto bottom = intersections[i]; if ( top.z() >= tvd && bottom.z() < tvd ) { topIdx = static_cast( i ) - 1; bottomIdx = static_cast( i ); break; } } } return std::make_pair( topIdx, bottomIdx ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::pair RimFaultReactivationDataAccessorWellLogExtraction::findOverburdenAndUnderburdenIndex( const std::vector& values ) { auto findLastOverburdenIndex = []( const std::vector& values ) { for ( size_t i = 0; i < values.size(); i++ ) { if ( !std::isinf( values[i] ) ) return static_cast( i ); } return -1; }; auto findFirstUnderburdenIndex = []( const std::vector& values ) { for ( size_t i = values.size() - 1; i > 0; i-- ) { if ( !std::isinf( values[i] ) ) return static_cast( i ); } return -1; }; int lastOverburdenIndex = findLastOverburdenIndex( values ); int firstUnderburdenIndex = findFirstUnderburdenIndex( values ); return { lastOverburdenIndex, firstUnderburdenIndex }; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- double RimFaultReactivationDataAccessorWellLogExtraction::computeValueWithGradient( const std::vector& intersections, const std::vector& values, int i1, int i2, double gradient ) { double tvdDiff = intersections[i2].z() - intersections[i1].z(); return tvdDiff * gradient + values[i2]; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValuesWithGradient( const std::vector& intersections, std::vector& values, double gradient ) { CAF_ASSERT( intersections.size() == values.size() ); auto calculatePorePressure = []( double depth, double gradient ) { return RiaEclipseUnitTools::pascalToBar( gradient * 9.81 * depth * 1000.0 ); }; auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values ); // Fill in overburden values using gradient double topPorePressure = calculatePorePressure( std::abs( intersections[0].z() ), gradient ); insertOverburdenValues( intersections, values, lastOverburdenIndex, topPorePressure ); // Fill in underburden values using gradient int lastElementIndex = static_cast( values.size() ) - 1; double bottomPorePressure = calculatePorePressure( std::abs( intersections[lastElementIndex].z() ), gradient ); insertUnderburdenValues( intersections, values, firstUnderburdenIndex, bottomPorePressure ); // Interpolate the missing values (should only be intra-reservoir by now) std::vector intersectionsZ = extractDepthValues( intersections ); RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValuesWithTopValue( const std::vector& intersections, std::vector& values, double topValue ) { CAF_ASSERT( intersections.size() == values.size() ); auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values ); // Fill in overburden values using gradient insertOverburdenValues( intersections, values, lastOverburdenIndex, topValue ); // Fill in underburden values double bottomValue = values[firstUnderburdenIndex - 1]; insertUnderburdenValues( intersections, values, firstUnderburdenIndex, bottomValue ); // Interpolate the missing values (should only be intra-reservoir by now) std::vector intersectionsZ = extractDepthValues( intersections ); RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, double seabedDepth, double bottomDepth, const cvf::Vec3d& offset ) { cvf::Vec3d faultTop = faultTopPosition + offset; cvf::Vec3d seabed( faultTop.x(), faultTop.y(), seabedDepth ); cvf::Vec3d faultBottom = faultBottomPosition + offset; cvf::Vec3d underburdenBottom( faultBottom.x(), faultBottom.y(), bottomDepth ); return { seabed, faultTop, faultBottom, underburdenBottom }; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RimFaultReactivationDataAccessorWellLogExtraction::generateMds( const std::vector& points ) { CAF_ASSERT( points.size() >= 2 ); // Assume first at zero, all other points relative to that. std::vector mds = { 0.0 }; double sum = 0.0; for ( size_t i = 1; i < points.size(); i++ ) { sum += points[i - 1].pointDistance( points[i] ); mds.push_back( sum ); } return mds; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::pair>, std::map>> RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( const RigFaultReactivationModel& model, RigEclipseCaseData& eclipseCaseData, double seabedDepth ) { auto [faultTopPosition, faultBottomPosition] = model.faultTopBottom(); auto faultNormal = model.faultNormal() ^ cvf::Vec3d::Z_AXIS; faultNormal.normalize(); double distanceFromFault = 1.0; auto [topDepth, bottomDepth] = model.depthTopBottom(); std::map> wellPaths; std::map> extractors; for ( auto gridPart : model.allGridParts() ) { double sign = model.normalPointsAt() == gridPart ? -1.0 : 1.0; std::vector wellPoints = RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, faultBottomPosition, seabedDepth, bottomDepth, sign * faultNormal * distanceFromFault ); cvf::ref wellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); wellPaths[gridPart] = wellPath; std::string errorName = "fault reactivation data access"; cvf::ref extractor = new RigEclipseWellLogExtractor( &eclipseCaseData, wellPath.p(), errorName ); extractors[gridPart] = extractor; } return { wellPaths, extractors }; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::pair, std::vector> RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( const RigResultAccessor& resultAccessor, RigEclipseWellLogExtractor& extractor, const RigWellPath& wellPath ) { // Extract values along well path std::vector values; extractor.curveData( &resultAccessor, &values ); auto intersections = extractor.intersections(); // Insert top of overburden point intersections.insert( intersections.begin(), wellPath.wellPathPoints().front() ); values.insert( values.begin(), std::numeric_limits::infinity() ); // Insert bottom of underburden point intersections.push_back( wellPath.wellPathPoints().back() ); values.push_back( std::numeric_limits::infinity() ); return { values, intersections }; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- double RimFaultReactivationDataAccessorWellLogExtraction::computeGradient( double depth1, double value1, double depth2, double value2 ) { return ( value2 - value1 ) / ( depth2 - depth1 ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RimFaultReactivationDataAccessorWellLogExtraction::extractDepthValues( const std::vector& intersections ) { std::vector intersectionsZ; for ( auto i : intersections ) { intersectionsZ.push_back( -i.z() ); } return intersectionsZ; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFaultReactivationDataAccessorWellLogExtraction::insertUnderburdenValues( const std::vector& intersections, std::vector& values, int firstUnderburdenIndex, double bottomValue ) { int lastElementIndex = static_cast( values.size() ) - 1; double underburdenGradient = computeGradient( intersections[firstUnderburdenIndex].z(), values[firstUnderburdenIndex], intersections[lastElementIndex].z(), bottomValue ); for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- ) { values[i] = computeValueWithGradient( intersections, values, i, firstUnderburdenIndex, -underburdenGradient ); } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFaultReactivationDataAccessorWellLogExtraction::insertOverburdenValues( const std::vector& intersections, std::vector& values, int lastOverburdenIndex, double topValue ) { double overburdenGradient = computeGradient( intersections[0].z(), topValue, intersections[lastOverburdenIndex].z(), values[lastOverburdenIndex] ); for ( int i = 0; i < lastOverburdenIndex; i++ ) { values[i] = computeValueWithGradient( intersections, values, i, lastOverburdenIndex, -overburdenGradient ); } }