Merge pull request #5199 from OPM/msw-fixes

Several critical fixes for MSW export

ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportMswCompletionsImpl.cpp

@@ -1043,40 +1043,12 @@ RicMswExportInfo RicWellPathExportMswCompletionsImpl::generatePerforationsMswExp
     int                                               timeStep,
     const std::vector<const RimPerforationInterval*>& perforationIntervals )
 {
-    const RigActiveCellInfo* activeCellInfo = eclipseCase->eclipseCaseData()->activeCellInfo( RiaDefines::MATRIX_MODEL );
     RiaEclipseUnitTools::UnitSystem unitSystem = eclipseCase->eclipseCaseData()->unitsType();
 
-    const RigWellPath*             wellPathGeometry = wellPath->wellPathGeometry();
+    double initialMD = 0.0; // Start measured depth location to export MSW data for. Either based on first intersection
-    const std::vector<cvf::Vec3d>& coords           = wellPathGeometry->wellPathPoints();
+                            // with active grid, or user defined value.
-    const std::vector<double>&     mds              = wellPathGeometry->measureDepths();
-    CVF_ASSERT( !coords.empty() && !mds.empty() );
 
-    std::vector<WellPathCellIntersectionInfo> intersections =
+    std::vector<SubSegmentIntersectionInfo> subSegIntersections = generateSubSegments( eclipseCase, wellPath, initialMD );
-        RigWellPathIntersectionTools::findCellIntersectionInfosAlongPath( eclipseCase->eclipseCaseData(), coords, mds );
-
-    double maxSegmentLength = wellPath->perforationIntervalCollection()->mswParameters()->maxSegmentLength();
-    std::vector<SubSegmentIntersectionInfo> subSegIntersections =
-        SubSegmentIntersectionInfo::spiltIntersectionSegmentsToMaxLength( wellPathGeometry,
-                                                                          intersections,
-                                                                          maxSegmentLength );
-
-    double initialMD = 0.0;
-    if ( wellPath->perforationIntervalCollection()->mswParameters()->referenceMDType() ==
-         RimMswCompletionParameters::MANUAL_REFERENCE_MD )
-    {
-        initialMD = wellPath->perforationIntervalCollection()->mswParameters()->manualReferenceMD();
-    }
-    else
-    {
-        for ( WellPathCellIntersectionInfo intersection : intersections )
-        {
-            if ( activeCellInfo->isActive( intersection.globCellIndex ) )
-            {
-                initialMD = intersection.startMD;
-                break;
-            }
-        }
-    }
 
     RicMswExportInfo exportInfo( wellPath,
                                  unitSystem,
@@ -1114,6 +1086,128 @@ RicMswExportInfo RicWellPathExportMswCompletionsImpl::generatePerforationsMswExp
     return exportInfo;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<SubSegmentIntersectionInfo>
+    RicWellPathExportMswCompletionsImpl::generateSubSegments( const RimEclipseCase* eclipseCase,
+                                                              const RimWellPath*    wellPath,
+                                                              double&               initialMD )
+{
+    const RigActiveCellInfo* activeCellInfo = eclipseCase->eclipseCaseData()->activeCellInfo( RiaDefines::MATRIX_MODEL );
+    const RigWellPath*       wellPathGeometry = wellPath->wellPathGeometry();
+    const std::vector<cvf::Vec3d>& coords     = wellPathGeometry->wellPathPoints();
+    const std::vector<double>&     mds        = wellPathGeometry->measureDepths();
+    CVF_ASSERT( !coords.empty() && !mds.empty() );
+
+    std::vector<WellPathCellIntersectionInfo> intersections =
+        RigWellPathIntersectionTools::findCellIntersectionInfosAlongPath( eclipseCase->eclipseCaseData(), coords, mds );
+
+    if ( wellPath->perforationIntervalCollection()->mswParameters()->referenceMDType() ==
+         RimMswCompletionParameters::MANUAL_REFERENCE_MD )
+    {
+        initialMD = wellPath->perforationIntervalCollection()->mswParameters()->manualReferenceMD();
+    }
+    else
+    {
+        for ( const WellPathCellIntersectionInfo& intersection : intersections )
+        {
+            if ( activeCellInfo->isActive( intersection.globCellIndex ) )
+            {
+                initialMD = intersection.startMD;
+                break;
+            }
+        }
+    }
+
+    std::vector<WellPathCellIntersectionInfo> filteredIntersections = filterIntersections( intersections,
+                                                                                           initialMD,
+                                                                                           wellPathGeometry,
+                                                                                           eclipseCase );
+
+    const double maxSegmentLength = wellPath->perforationIntervalCollection()->mswParameters()->maxSegmentLength();
+
+    return SubSegmentIntersectionInfo::spiltIntersectionSegmentsToMaxLength( wellPathGeometry,
+                                                                             filteredIntersections,
+                                                                             maxSegmentLength );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<WellPathCellIntersectionInfo> RicWellPathExportMswCompletionsImpl::filterIntersections(
+    const std::vector<WellPathCellIntersectionInfo>& intersections,
+    double                                           initialMD,
+    const RigWellPath*                               wellPathGeometry,
+    const RimEclipseCase*                            eclipseCase )
+{
+    std::vector<WellPathCellIntersectionInfo> filteredIntersections;
+
+    if ( !intersections.empty() && intersections[0].startMD > initialMD )
+    {
+        WellPathCellIntersectionInfo firstIntersection = intersections[0];
+
+        // Add a segment from user defined MD to start of grid
+        cvf::Vec3d intersectionPoint = wellPathGeometry->interpolatedPointAlongWellPath( initialMD );
+
+        WellPathCellIntersectionInfo extraIntersection;
+
+        extraIntersection.globCellIndex          = std::numeric_limits<size_t>::max();
+        extraIntersection.startPoint             = intersectionPoint;
+        extraIntersection.endPoint               = firstIntersection.startPoint;
+        extraIntersection.startMD                = initialMD;
+        extraIntersection.endMD                  = firstIntersection.startMD;
+        extraIntersection.intersectedCellFaceIn  = cvf::StructGridInterface::NO_FACE;
+        extraIntersection.intersectedCellFaceOut = cvf::StructGridInterface::oppositeFace(
+            firstIntersection.intersectedCellFaceIn );
+        extraIntersection.intersectionLengthsInCellCS = cvf::Vec3d::ZERO;
+
+        filteredIntersections.push_back( extraIntersection );
+    }
+
+    const double epsilon = 1.0e-3;
+
+    for ( const WellPathCellIntersectionInfo& intersection : intersections )
+    {
+        if ( ( intersection.endMD - initialMD ) < epsilon )
+        {
+            // Skip all intersections before initial measured depth
+            continue;
+        }
+        else if ( ( intersection.startMD - initialMD ) > epsilon )
+        {
+            filteredIntersections.push_back( intersection );
+        }
+        else
+        {
+            // InitialMD is inside intersection, split based on intersection point
+
+            cvf::Vec3d intersectionPoint = wellPathGeometry->interpolatedPointAlongWellPath( initialMD );
+
+            WellPathCellIntersectionInfo extraIntersection;
+
+            extraIntersection.globCellIndex          = intersection.globCellIndex;
+            extraIntersection.startPoint             = intersectionPoint;
+            extraIntersection.endPoint               = intersection.endPoint;
+            extraIntersection.startMD                = initialMD;
+            extraIntersection.endMD                  = intersection.endMD;
+            extraIntersection.intersectedCellFaceIn  = cvf::StructGridInterface::NO_FACE;
+            extraIntersection.intersectedCellFaceOut = intersection.intersectedCellFaceOut;
+
+            const RigMainGrid* grid = eclipseCase->mainGrid();
+
+            extraIntersection.intersectionLengthsInCellCS =
+                RigWellPathIntersectionTools::calculateLengthInCell( grid,
+                                                                     intersection.globCellIndex,
+                                                                     intersectionPoint,
+                                                                     intersection.endPoint );
+            filteredIntersections.push_back( extraIntersection );
+        }
+    }
+
+    return filteredIntersections;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -1128,9 +1222,16 @@ RicWellPathExportMswCompletionsImpl::MainBoreSegments
 {
     MainBoreSegments mainBoreSegments;
 
+    // Intersections along the well path with grid geometry is handled by well log extraction tools. The threshold in
+    // RigWellLogExtractionTools::isEqualDepth is currently set to 0.1m, and this is a pretty large threshold based on
+    // the indicated threshold of 0.001m for MSW segments
+    const double segmentLengthThreshold = 1.0e-3;
+
     for ( const auto& cellIntInfo : subSegIntersections )
     {
-        if ( std::fabs( cellIntInfo.endMD - cellIntInfo.startMD ) > 1.0e-8 )
+        const double segmentLength = std::fabs( cellIntInfo.endMD - cellIntInfo.startMD );
+
+        if ( segmentLength > segmentLengthThreshold )
         {
             QString                        label = QString( "Main stem segment %1" ).arg( mainBoreSegments.size() + 2 );
             std::shared_ptr<RicMswSegment> segment( new RicMswSegment( label,
@@ -1162,6 +1263,13 @@ RicWellPathExportMswCompletionsImpl::MainBoreSegments
             }
             mainBoreSegments.push_back( segment );
         }
+        else
+        {
+            QString text = QString( "Skipping segment , threshold = %1, length = %2" )
+                               .arg( segmentLengthThreshold )
+                               .arg( segmentLength );
+            RiaLogging::info( text );
+        }
     }
     return mainBoreSegments;
 }

ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportMswCompletionsImpl.h

@@ -29,6 +29,9 @@ class RimWellPath;
 class RimWellPathValve;
 class RimWellPathFracture;
 class SubSegmentIntersectionInfo;
+class RigWellPath;
+
+struct WellPathCellIntersectionInfo;
 
 class QFile;
 
@@ -76,6 +79,15 @@ private:
                                            int                                               timeStep,
                                            const std::vector<const RimPerforationInterval*>& perforationIntervals );
 
+    static std::vector<SubSegmentIntersectionInfo>
+        generateSubSegments( const RimEclipseCase* eclipseCase, const RimWellPath* wellPath, double& initialMD );
+
+    static std::vector<WellPathCellIntersectionInfo>
+        filterIntersections( const std::vector<WellPathCellIntersectionInfo>& intersections,
+                             double                                           initialMD,
+                             const RigWellPath*                               wellPathGeometry,
+                             const RimEclipseCase*                            eclipseCase );
+
     static void generateWelsegsTable( RifTextDataTableFormatter& formatter, const RicMswExportInfo& exportInfo );
 
     static void generateWelsegsSegments( RifTextDataTableFormatter&                         formatter,