feat: Add functions to connect well perforation to surface elements

inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml → 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>

inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_smoke.xml → 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
@@ -57,15 +58,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>
 

inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml → 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
@@ -39,6 +40,7 @@
       targetRegions="{ Region, Fault, wellRegion1, wellRegion2 }">
       <NonlinearSolverParameters
         newtonMaxIter="40"
+        newtonTol="1e-4"
         maxAllowedResidualNorm="1e+25"/>
       <LinearSolverParameters
         directParallel="0"/>
@@ -61,15 +63,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>
 

src/coreComponents/mainInterface/ProblemManager.cpp

@@ -865,7 +865,7 @@ void ProblemManager::generateMeshLevel( MeshLevel & meshLevel,
     // 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 ) )
+    if( isBaseMeshLevel ) // && !dynamicCast< FaceElementSubRegion * >( &subRegion ) )
     {
       subRegion.calculateElementGeometricQuantities( nodeManager, faceManager );
     }

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 )) )

src/coreComponents/mesh/ElementRegionManager.cpp

@@ -204,7 +204,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();

src/coreComponents/mesh/ElementSubRegionBase.hpp

@@ -280,13 +280,23 @@ class ElementSubRegionBase : public ObjectManagerBase
     forAll< parallelHostPolicy >( size(), [=]( localIndex const k )
     {
       LvArray::tensorOps::copy< 3 >( elementCenters[k], X[e2n( k, 0 )] );
+      std::set< std::array< real64, 3 > > uniqueNodes;
+      {
+        std::array< real64, 3 > const node = {X( e2n( k, 0 ), 0 ), X( e2n( k, 0 ), 1 ), X( e2n( k, 0 ), 2 )};
+        uniqueNodes.insert( node );
+      }
       localIndex const numNodes = this->numNodesPerElement( k );
       for( localIndex a = 1; a < numNodes; ++a )
       {
-        LvArray::tensorOps::add< 3 >( elementCenters[k], X[e2n( k, a )] );
+        std::array< real64, 3 > const node = {X[e2n( k, a )][0], X[e2n( k, a )][1], X[e2n( k, a )][2]};
+        if( uniqueNodes.find( node ) == uniqueNodes.end())
+        {
+          uniqueNodes.insert( node );
+          LvArray::tensorOps::add< 3 >( elementCenters[k], X[e2n( k, a )] );
+        }
       }
 
-      LvArray::tensorOps::scale< 3 >( elementCenters[k], 1.0 / numNodes );
+      LvArray::tensorOps::scale< 3 >( elementCenters[k], 1.0 / uniqueNodes.size() );
     } );
   }
 };

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." );

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

src/coreComponents/mesh/PerforationData.cpp

@@ -233,23 +233,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,

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

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 );
 
   ///@}