feat: Add functions to connect well perforation to surface elements (#3359)

paveltomin · jhuang2601 · web-flow · commit 0f7eab56e4c5 · 2025-10-07T13:49:20.000-05:00
Co-authored-by: Jian HUANG &lt;jian.huang@totalenergies.com&gt;
diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3843-13975-76d125d
+  baseline: integratedTests/baseline_integratedTests-pr3359-14024-3c9ebb4
 
 allow_fail:
   all: ''
diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md
@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3359 (2025-10-07) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3359-14024-3c9ebb4.tar.gz>
+=====================
+Add functions to connect well perforation to surface elements.
+
 PR #3843 (2025-10-02) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3843-13975-76d125d.tar.gz>
 =====================
 Fix linear solver issues.
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_base.xml
@@ -29,7 +29,7 @@
 
     <ConstantPermeability
       name="rockPerm"
-      permeabilityComponents="{ 1.0e-18, 1.0e-18, 1.0e-18 }" />
+      permeabilityComponents="{ 1.0e-14, 1.0e-14, 1.0e-14 }" />
       <!-- Material inside the fault -->
     <CompressibleSolidParallelPlatesPermeability
       name="faultFilling"
@@ -188,14 +188,15 @@
         logLevel="1"
         wellRegionName="wellRegion2"
         wellControlsName="wellControls2"
-        polylineNodeCoords="{ { -5400, -5400, 0 },
-                              { -5400, -5400, -2500 } }"
+        polylineNodeCoords="{ { 1200, -700, 0 },
+                              { 1200, -700, -2500 } }"
         polylineSegmentConn="{ { 0, 1 } }"
         radius="0.1"
         numElementsPerSegment="1">
         <Perforation
           name="injector1_perf1"
-          distanceFromHead="2500"/>
+          targetRegion="Fault"
+          distanceFromHead="1700"/>
       </InternalWell>
     </VTKMesh>
   </Mesh>
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_fim_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_fim_smoke.xml
@@ -1,6 +1,6 @@
 <Problem>
   <Included>
-    <File name="./singlePhasePoromechanics_FaultModel_well_base.xml"/>
+    <File name="./singlePhasePoromechanics_WellInFault_base.xml"/>
   </Included>
 
   <Solvers
@@ -20,7 +20,8 @@
         newtonMaxIter="20"
         maxAllowedResidualNorm="1e+25"/>
       <LinearSolverParameters
-        solverType="direct"/> 
+        solverType="direct"
+        directParallel="1"/> 
     </SinglePhasePoromechanicsConformingFracturesReservoir>
 
     <SinglePhasePoromechanicsConformingFractures
@@ -56,15 +57,15 @@
         control="BHP"
         referenceElevation="-2500"
         targetBHP="1e7"
-        targetTotalRate="1e10" />
+        targetTotalRate="1e-3" />
       <WellControls
         name="wellControls2"
         logLevel="2"
         type="injector"
         control="BHP"
         referenceElevation="-2500"
         targetBHP="10e7"
-        targetTotalRate="1e10" />
+        targetTotalRate="1e-3" />
     </SinglePhaseWell>
   </Solvers>
 
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_seq_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_seq_smoke.xml
@@ -1,6 +1,6 @@
 <Problem>
   <Included>
-    <File name="./singlePhasePoromechanics_FaultModel_well_new_base.xml"/>
+    <File name="./singlePhasePoromechanics_WellInFault_base.xml"/>
   </Included>
 
   <Solvers
@@ -20,7 +20,8 @@
         newtonTol="1e-3"
         newtonMaxIter="20"/>
       <LinearSolverParameters
-        solverType="direct"/> 
+        solverType="direct"
+        directParallel="1"/> 
     </SinglePhaseReservoirPoromechanicsConformingFractures>
 
     <SolidMechanicsLagrangeContact
@@ -38,6 +39,7 @@
       targetRegions="{ Region, Fault, wellRegion1, wellRegion2 }">
       <NonlinearSolverParameters
         newtonMaxIter="40"
+        newtonTol="1e-4"
         maxAllowedResidualNorm="1e+25"/>
       <LinearSolverParameters
         directParallel="0"/>
@@ -60,15 +62,15 @@
         control="BHP"
         referenceElevation="-2500"
         targetBHP="1e7"
-        targetTotalRate="1e10" />
+        targetTotalRate="1e-3" />
       <WellControls
         name="wellControls2"
         logLevel="2"
         type="injector"
         control="BHP"
         referenceElevation="-2500"
         targetBHP="10e7"
-        targetTotalRate="1e10" />
+        targetTotalRate="1e-3" />
     </SinglePhaseWell>
   </Solvers>
 
diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp
@@ -863,11 +863,12 @@ void ProblemManager::generateMeshLevel( MeshLevel & meshLevel,
   elemRegionManager.forElementSubRegions< ElementSubRegionBase >( [&]( ElementSubRegionBase & subRegion )
   {
     subRegion.setupRelatedObjectsInRelations( meshLevel );
+    // TODO calling calculateElementGeometricQuantities for `FaceElementSubRegion` here is not very accurate:
     // `FaceElementSubRegion` has no node and therefore needs the nodes positions from the neighbor elements
     // in order to compute the geometric quantities.
     // And this point of the process, the ghosting has not been done and some elements of the `FaceElementSubRegion`
-    // can have no neighbor. Making impossible the computation, which is therfore postponed to after the ghosting.
-    if( isBaseMeshLevel && !dynamicCast< FaceElementSubRegion * >( &subRegion ) )
+    // can have no neighbor. We still call it only to compute element centers to be used for well perforation computations.
+    if( isBaseMeshLevel ) // && !dynamicCast< FaceElementSubRegion * >( &subRegion ) )
     {
       subRegion.calculateElementGeometricQuantities( nodeManager, faceManager );
     }
diff --git a/src/coreComponents/mesh/DomainPartition.cpp b/src/coreComponents/mesh/DomainPartition.cpp
@@ -247,9 +247,10 @@ void DomainPartition::setupCommunications( bool use_nonblocking )
         {
           NodeManager & nodeManager = meshLevel.getNodeManager();
           FaceManager & faceManager = meshLevel.getFaceManager();
+          ElementRegionManager & elemManager = meshLevel.getElemManager();
 
           CommunicationTools::getInstance().setupGhosts( meshLevel, m_neighbors, use_nonblocking );
-          faceManager.sortAllFaceNodes( nodeManager, meshLevel.getElemManager() );
+          faceManager.sortAllFaceNodes( nodeManager, elemManager );
           faceManager.computeGeometry( nodeManager );
         }
         else if( !meshLevel.isShallowCopyOf( meshBody.getMeshLevels().getGroup< MeshLevel >( 0 )) )
diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp
@@ -203,7 +203,11 @@ void ElementRegionManager::generateWells( CellBlockManagerABC const & cellBlockM
     // generate the local data (well elements, nodes, perforations) on this well
     // note: each MPI rank knows the global info on the entire well (constructed earlier in InternalWellGenerator)
     // so we only need node and element offsets to construct the local-to-global maps in each wellElemSubRegion
-    wellRegion.generateWell( meshLevel, lineBlock, nodeOffsetGlobal + wellNodeCount, elemOffsetGlobal + wellElemCount );
+    wellRegion.generateWell( meshLevel,
+                             lineBlock,
+                             nodeOffsetGlobal + wellNodeCount,
+                             elemOffsetGlobal + wellElemCount,
+                             cellBlockManager.getGlobalLength() );
 
     // increment counters with global number of nodes and elements
     wellElemCount += lineBlock.numElements();
diff --git a/src/coreComponents/mesh/ElementSubRegionBase.hpp b/src/coreComponents/mesh/ElementSubRegionBase.hpp
@@ -279,14 +279,45 @@ class ElementSubRegionBase : public ObjectManagerBase
 
     forAll< parallelHostPolicy >( size(), [=]( localIndex const k )
     {
-      LvArray::tensorOps::copy< 3 >( elementCenters[k], X[e2n( k, 0 )] );
+      // collect node coordinates for element k
       localIndex const numNodes = this->numNodesPerElement( k );
-      for( localIndex a = 1; a < numNodes; ++a )
+      std::vector< std::array< real64, 3 > > nodes;
+      nodes.reserve( numNodes );
+      for( localIndex a = 0; a < numNodes; ++a )
       {
-        LvArray::tensorOps::add< 3 >( elementCenters[k], X[e2n( k, a )] );
+        localIndex const nA = e2n( k, a );
+        nodes.push_back( { X( nA, 0 ), X( nA, 1 ), X( nA, 2 ) } );
       }
 
-      LvArray::tensorOps::scale< 3 >( elementCenters[k], 1.0 / numNodes );
+      // sort + unique with tolerance to remove near-duplicates (stable and robust for floating point)
+      auto cmpLex = []( auto const & A, auto const & B )
+      {
+        return std::tie( A[0], A[1], A[2] ) < std::tie( B[0], B[1], B[2] );
+      };
+      const real64 tol = 1e-10;
+      auto almostEqual = [tol]( auto const & A, auto const & B )
+      {
+        return ( std::abs( A[0] - B[0] ) < tol ) &&
+        ( std::abs( A[1] - B[1] ) < tol ) &&
+        ( std::abs( A[2] - B[2] ) < tol );
+      };
+
+      std::sort( nodes.begin(), nodes.end(), cmpLex );
+      nodes.erase( std::unique( nodes.begin(), nodes.end(), almostEqual ), nodes.end() );
+
+      // compute average (center) of unique nodes
+      std::array< real64, 3 > centre = { 0.0, 0.0, 0.0 };
+      for( auto const & p : nodes )
+      {
+        centre[0] += p[0];
+        centre[1] += p[1];
+        centre[2] += p[2];
+      }
+      const real64 inv = 1.0 / static_cast< real64 >( nodes.size() );
+
+      elementCenters[k][0] = centre[0] * inv;
+      elementCenters[k][1] = centre[1] * inv;
+      elementCenters[k][2] = centre[2] * inv;
     } );
   }
 };
diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp
@@ -48,14 +48,6 @@ EmbeddedSurfaceSubRegion::EmbeddedSurfaceSubRegion( string const & name,
 {
   m_elementType = ElementType::Polygon;
 
-  registerWrapper( viewKeyStruct::elementCenterString(), &m_elementCenter ).
-    setDescription( "The center of each EmbeddedSurface element." );
-
-  registerWrapper( viewKeyStruct::elementVolumeString(), &m_elementVolume ).
-    setApplyDefaultValue( -1.0 ).
-    setPlotLevel( dataRepository::PlotLevel::LEVEL_0 ).
-    setDescription( "The volume of each EmbeddedSurface element." );
-
   registerWrapper( viewKeyStruct::connectivityIndexString(), &m_connectivityIndex ).
     setApplyDefaultValue( 1 ).
     setDescription( "Connectivity index of each EmbeddedSurface." );
diff --git a/src/coreComponents/mesh/FaceElementSubRegion.cpp b/src/coreComponents/mesh/FaceElementSubRegion.cpp
@@ -230,7 +230,7 @@ void FaceElementSubRegion::calculateSingleElementGeometricQuantities( localIndex
   m_elementVolume[k] = m_elementAperture[k] * faceArea[m_toFacesRelation[k][0]];
 }
 
-void FaceElementSubRegion::calculateElementGeometricQuantities( NodeManager const & GEOS_UNUSED_PARAM( nodeManager ),
+void FaceElementSubRegion::calculateElementGeometricQuantities( NodeManager const & nodeManager,
                                                                 FaceManager const & faceManager )
 {
   arrayView1d< real64 const > const & faceArea = faceManager.faceArea();
@@ -239,6 +239,9 @@ void FaceElementSubRegion::calculateElementGeometricQuantities( NodeManager cons
   {
     calculateSingleElementGeometricQuantities( k, faceArea );
   } );
+
+  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X = nodeManager.referencePosition();
+  calculateElementCenters( X );
 }
 
 ElementType FaceElementSubRegion::getElementType( localIndex ei ) const
diff --git a/src/coreComponents/mesh/PerforationData.cpp b/src/coreComponents/mesh/PerforationData.cpp
@@ -246,23 +246,40 @@ void PerforationData::getReservoirElementDimensions( MeshLevel const & mesh,
 {
   ElementRegionManager const & elemManager = mesh.getElemManager();
   NodeManager const & nodeManager = mesh.getNodeManager();
-  CellElementRegion const & region = elemManager.getRegion< CellElementRegion >( er );
-  CellElementSubRegion const & subRegion = region.getSubRegion< CellElementSubRegion >( esr );
-
-  // compute the bounding box of the element
-  real64 boxDims[ 3 ];
-  computationalGeometry::getBoundingBox( ei,
-                                         subRegion.nodeList(),
-                                         nodeManager.referencePosition(),
-                                         boxDims );
-
-  // dx and dz from bounding box
-  dx = boxDims[ 0 ];
-  dy = boxDims[ 1 ];
-
-  // dz is computed as vol / (dx * dy)
-  dz  = subRegion.getElementVolume()[ei];
-  dz /= dx * dy;
+  ElementRegionBase const & region = elemManager.getRegion< ElementRegionBase >( er );
+  ElementSubRegionBase const & subRegionBase = region.getSubRegion< ElementSubRegionBase >( esr );
+
+  region.applyLambdaToContainer< CellElementSubRegion, SurfaceElementSubRegion >( subRegionBase, [&]( auto const & subRegion )
+  {
+    // compute the bounding box of the element
+    real64 boxDims[ 3 ];
+    computationalGeometry::getBoundingBox( ei,
+                                           subRegion.nodeList(),
+                                           nodeManager.referencePosition(),
+                                           boxDims );
+
+    // dx and dz from bounding box
+    dx = boxDims[ 0 ];
+    dy = boxDims[ 1 ];
+    dz = boxDims[ 2 ];
+
+    if( isZero( dx ) )
+    {
+      dx  = subRegion.getElementVolume()[ei];
+      dx /= dy * dz;
+    }
+    else if( isZero( dy ) )
+    {
+      dy  = subRegion.getElementVolume()[ei];
+      dy /= dx * dz;
+    }
+    else
+    {
+      // dz is computed as vol / (dx * dy)
+      dz  = subRegion.getElementVolume()[ei];
+      dz /= dx * dy;
+    }
+  } );
 }
 
 void PerforationData::connectToWellElements( LineBlockABC const & lineBlock,
diff --git a/src/coreComponents/mesh/WellElementRegion.cpp b/src/coreComponents/mesh/WellElementRegion.cpp
@@ -49,7 +49,8 @@ WellElementRegion::~WellElementRegion()
 void WellElementRegion::generateWell( MeshLevel & mesh,
                                       LineBlockABC const & lineBlock,
                                       globalIndex nodeOffsetGlobal,
-                                      globalIndex elemOffsetGlobal )
+                                      globalIndex elemOffsetGlobal,
+                                      real64 const globalLength )
 {
   // get the (unique) subregion
   WellElementSubRegion &
@@ -64,8 +65,11 @@ void WellElementRegion::generateWell( MeshLevel & mesh,
   globalIndex const numElemsGlobal        = lineBlock.numElements();
   globalIndex const numPerforationsGlobal = lineBlock.numPerforations();
 
+  // tolerance for geometrical calculations
+  real64 const geomTol = globalLength * 1e-10;
+
   // 1) select the local perforations based on connectivity to the local reservoir elements
-  subRegion.connectPerforationsToMeshElements( mesh, lineBlock );
+  subRegion.connectPerforationsToMeshElements( mesh, lineBlock, geomTol );
 
   globalIndex const matchedPerforations = MpiWrapper::sum( perforationData->size() );
   GEOS_THROW_IF( matchedPerforations != numPerforationsGlobal,
@@ -105,7 +109,8 @@ void WellElementRegion::generateWell( MeshLevel & mesh,
                       lineBlock,
                       elemStatusGlobal,
                       nodeOffsetGlobal,
-                      elemOffsetGlobal );
+                      elemOffsetGlobal,
+                      geomTol );
 
 
   // 4) find out which rank is the owner of the top segment
diff --git a/src/coreComponents/mesh/WellElementRegion.hpp b/src/coreComponents/mesh/WellElementRegion.hpp
@@ -130,11 +130,13 @@ class WellElementRegion : public ElementRegionBase
    * @param[in] lineBlock the LineBlockABC containing the global well topology
    * @param[in] nodeOffsetGlobal the offset of the first global well node ( = offset of last global mesh node + 1 )
    * @param[in] elemOffsetGlobal the offset of the first global well element ( = offset of last global mesh elem + 1 )
+   * @param[in] globalLength the global length of the mesh
    */
   void generateWell( MeshLevel & mesh,
                      LineBlockABC const & lineBlock,
                      globalIndex nodeOffsetGlobal,
-                     globalIndex elemOffsetGlobal );
+                     globalIndex elemOffsetGlobal,
+                     real64 globalLength );
 
   ///@}
 
diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp
diff --git a/src/coreComponents/mesh/WellElementSubRegion.hpp b/src/coreComponents/mesh/WellElementSubRegion.hpp
diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp

Original file line number	Diff line number	Diff line change
`@@ -863,11 +863,12 @@ void ProblemManager::generateMeshLevel( MeshLevel & meshLevel,`
`863`	`863`	`elemRegionManager.forElementSubRegions< ElementSubRegionBase >( [&]( ElementSubRegionBase & subRegion )`
`864`	`864`	`{`
`865`	`865`	`subRegion.setupRelatedObjectsInRelations( meshLevel );`
	`866`	+ // TODO calling calculateElementGeometricQuantities for `FaceElementSubRegion` here is not very accurate:
`866`	`867`	// `FaceElementSubRegion` has no node and therefore needs the nodes positions from the neighbor elements
`867`	`868`	`// in order to compute the geometric quantities.`
`868`	`869`	// And this point of the process, the ghosting has not been done and some elements of the `FaceElementSubRegion`
`869`		`- // can have no neighbor. Making impossible the computation, which is therfore postponed to after the ghosting.`
`870`		`- if( isBaseMeshLevel && !dynamicCast< FaceElementSubRegion * >( &subRegion ) )`
	`870`	`+ // can have no neighbor. We still call it only to compute element centers to be used for well perforation computations.`
	`871`	`+ if( isBaseMeshLevel ) // && !dynamicCast< FaceElementSubRegion * >( &subRegion ) )`
`871`	`872`	`{`
`872`	`873`	`subRegion.calculateElementGeometricQuantities( nodeManager, faceManager );`
`873`	`874`	`}`
Original file line number	Diff line number	Diff line change
`@@ -247,9 +247,10 @@ void DomainPartition::setupCommunications( bool use_nonblocking )`
`247`	`247`	`{`
`248`	`248`	`NodeManager & nodeManager = meshLevel.getNodeManager();`
`249`	`249`	`FaceManager & faceManager = meshLevel.getFaceManager();`
	`250`	`+ ElementRegionManager & elemManager = meshLevel.getElemManager();`
`250`	`251`
`251`	`252`	`CommunicationTools::getInstance().setupGhosts( meshLevel, m_neighbors, use_nonblocking );`
`252`		`- faceManager.sortAllFaceNodes( nodeManager, meshLevel.getElemManager() );`
	`253`	`+ faceManager.sortAllFaceNodes( nodeManager, elemManager );`
`253`	`254`	`faceManager.computeGeometry( nodeManager );`
`254`	`255`	`}`
`255`	`256`	`else if( !meshLevel.isShallowCopyOf( meshBody.getMeshLevels().getGroup< MeshLevel >( 0 )) )`