Fixed source injection when using multi region with first order wave …

…propagation solver (GEOS-DEV#2488)

* Fixed bug with sources injection + receivers storage when we use several regions

src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp

@@ -70,6 +70,16 @@ AcousticFirstOrderWaveEquationSEM::AcousticFirstOrderWaveEquationSEM( const std:
     setSizedFromParent( 0 ).
     setDescription( "Element containing the receivers" );
 
+  registerWrapper( viewKeyStruct::sourceRegionString(), &m_sourceRegion ).
+    setInputFlag( InputFlags::FALSE ).
+    setSizedFromParent( 0 ).
+    setDescription( "Region containing the sources" );
+
+  registerWrapper( viewKeyStruct::receiverRegionString(), &m_receiverRegion ).
+    setInputFlag( InputFlags::FALSE ).
+    setSizedFromParent( 0 ).
+    setDescription( "Region containing the receivers" );
+
 }
 
 AcousticFirstOrderWaveEquationSEM::~AcousticFirstOrderWaveEquationSEM()
@@ -145,9 +155,9 @@ void AcousticFirstOrderWaveEquationSEM::postProcessInput()
   m_uyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
   m_uzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
   m_sourceElem.resize( numSourcesGlobal );
+  m_sourceRegion.resize( numSourcesGlobal );
   m_rcvElem.resize( numReceiversGlobal );
-
-
+  m_receiverRegion.resize( numReceiversGlobal );
 }
 
 void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames )
@@ -165,6 +175,7 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
   arrayView2d< real64 > const sourceConstants = m_sourceConstants.toView();
   arrayView1d< localIndex > const sourceIsAccessible = m_sourceIsAccessible.toView();
   arrayView1d< localIndex > const sourceElem = m_sourceElem.toView();
+  arrayView1d< localIndex > const sourceRegion = m_sourceRegion.toView();
   sourceNodeIds.setValues< EXEC_POLICY >( -1 );
   sourceConstants.setValues< EXEC_POLICY >( -1 );
   sourceIsAccessible.zero();
@@ -174,6 +185,7 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
   arrayView2d< real64 > const receiverConstants = m_receiverConstants.toView();
   arrayView1d< localIndex > const receiverIsLocal = m_receiverIsLocal.toView();
   arrayView1d< localIndex > const rcvElem = m_rcvElem.toView();
+  arrayView1d< localIndex > const receiverRegion = m_receiverRegion.toView();
   receiverNodeIds.setValues< EXEC_POLICY >( -1 );
   receiverConstants.setValues< EXEC_POLICY >( -1 );
   receiverIsLocal.zero();
@@ -192,8 +204,7 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
     }
   }
 
-
-  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
+  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const regionIndex,
                                                                                         CellElementSubRegion & elementSubRegion )
   {
     GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron,
@@ -218,6 +229,7 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
         PrecomputeSourceAndReceiverKernel::
         launch< EXEC_POLICY, FE_TYPE >
         ( elementSubRegion.size(),
+        regionIndex,
         numNodesPerElem,
         numFacesPerElem,
         X,
@@ -232,11 +244,13 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
         sourceElem,
         sourceNodeIds,
         sourceConstants,
+        sourceRegion,
         receiverCoordinates,
         receiverIsLocal,
         rcvElem,
         receiverNodeIds,
         receiverConstants,
+        receiverRegion,
         sourceValue,
         dt,
         timeSourceFrequency,
@@ -439,6 +453,7 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
   arrayView2d< real64 const > const sourceConstants = m_sourceConstants.toView();
   arrayView1d< localIndex const > const sourceIsAccessible = m_sourceIsAccessible.toView();
   arrayView1d< localIndex const > const sourceElem = m_sourceElem.toView();
+  arrayView1d< localIndex const > const sourceRegion = m_sourceRegion.toView();
   arrayView2d< real32 const > const sourceValue = m_sourceValue.toView();
 
   GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc );
@@ -459,7 +474,7 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
 
     arrayView1d< real32 > const rhs = nodeManager.getField< wavesolverfields::ForcingRHS >();
 
-    mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
+    mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const regionIndex,
                                                                                           CellElementSubRegion & elementSubRegion )
     {
       arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList();
@@ -491,6 +506,7 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
           PressureComputation< FE_TYPE > kernel2( finiteElement );
         kernel2.template launch< EXEC_POLICY, ATOMIC_POLICY >
           ( elementSubRegion.size(),
+          regionIndex,
           nodeManager.size(),
           X,
           elemsToNodes,
@@ -503,6 +519,7 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
           sourceValue,
           sourceIsAccessible,
           sourceElem,
+          sourceRegion,
           dt,
           cycleNumber,
           p_np1 );
@@ -511,9 +528,9 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
       arrayView2d< real32 > const uyReceivers   = m_uyNp1AtReceivers.toView();
       arrayView2d< real32 > const uzReceivers   = m_uzNp1AtReceivers.toView();
 
-      compute2dVariableAllSeismoTraces( time_n, dt, velocity_x, velocity_x, uxReceivers );
-      compute2dVariableAllSeismoTraces( time_n, dt, velocity_y, velocity_y, uyReceivers );
-      compute2dVariableAllSeismoTraces( time_n, dt, velocity_z, velocity_z, uzReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_x, velocity_x, uxReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_y, velocity_y, uyReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_z, velocity_z, uzReceivers );
 
     } );
 
@@ -550,7 +567,7 @@ void AcousticFirstOrderWaveEquationSEM::cleanup( real64 const time_n, integer co
   {
     NodeManager & nodeManager = mesh.getNodeManager();
     arrayView1d< real32 const > const p_np1 = nodeManager.getField< wavesolverfields::Pressure_np1 >();
-    mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
+    mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const regionIndex,
                                                                                           CellElementSubRegion & elementSubRegion )
     {
       arrayView2d< real32 > const velocity_x = elementSubRegion.getField< wavesolverfields::Velocity_x >();
@@ -561,9 +578,9 @@ void AcousticFirstOrderWaveEquationSEM::cleanup( real64 const time_n, integer co
       arrayView2d< real32 > const uyReceivers   = m_uyNp1AtReceivers.toView();
       arrayView2d< real32 > const uzReceivers   = m_uzNp1AtReceivers.toView();
 
-      compute2dVariableAllSeismoTraces( time_n, 0, velocity_x, velocity_x, uxReceivers );
-      compute2dVariableAllSeismoTraces( time_n, 0, velocity_y, velocity_y, uyReceivers );
-      compute2dVariableAllSeismoTraces( time_n, 0, velocity_z, velocity_z, uzReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_x, velocity_x, uxReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_y, velocity_y, uyReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_z, velocity_z, uzReceivers );
 
     } );
     arrayView2d< real32 > const pReceivers   = m_pressureNp1AtReceivers.toView();
@@ -594,7 +611,8 @@ void AcousticFirstOrderWaveEquationSEM::computeAllSeismoTraces( real64 const tim
   }
 }
 
-void AcousticFirstOrderWaveEquationSEM::compute2dVariableAllSeismoTraces( real64 const time_n,
+void AcousticFirstOrderWaveEquationSEM::compute2dVariableAllSeismoTraces( localIndex const regionIndex,
+                                                                          real64 const time_n,
                                                                           real64 const dt,
                                                                           arrayView2d< real32 const > const var_np1,
                                                                           arrayView2d< real32 const > const var_n,
@@ -605,7 +623,7 @@ void AcousticFirstOrderWaveEquationSEM::compute2dVariableAllSeismoTraces( real64
        (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
        indexSeismoTrace++ )
   {
-    WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_rcvElem, m_receiverConstants, m_receiverIsLocal, m_nsamplesSeismoTrace,
+    WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dt, regionIndex, m_receiverRegion, timeSeismo, indexSeismoTrace, m_rcvElem, m_receiverConstants, m_receiverIsLocal, m_nsamplesSeismoTrace,
                                                    m_outputSeismoTrace,
                                                    var_np1, var_n, varAtReceivers );
   }

src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp

@@ -102,7 +102,8 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
    * @param var_n the field values at time_n
    * @param var_receivers the array holding the trace values, where the output is written
    */
-  virtual void compute2dVariableAllSeismoTraces( real64 const time_n,
+  virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex,
+                                                 real64 const time_n,
                                                  real64 const dt,
                                                  arrayView2d< real32 const > const var_np1,
                                                  arrayView2d< real32 const > const var_n,
@@ -131,7 +132,9 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
     static constexpr char const * uzNp1AtReceiversString() { return "uzNp1AtReceivers"; }
 
     static constexpr char const * sourceElemString() { return "sourceElem"; }
+    static constexpr char const * sourceRegionString() { return "sourceRegion"; }
     static constexpr char const * receiverElemString() { return "rcvElem"; }
+    static constexpr char const * receiverRegionString() { return "receiverRegion"; }
 
   } waveEquationViewKeys;
 
@@ -192,6 +195,12 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
   /// Array containing the elements which contain a source
   array1d< localIndex > m_sourceElem;
 
+  /// Array containing the elements which contain the region which the source belongs
+  array1d< localIndex > m_sourceRegion;
+
+  /// Array containing the elements which contain the region which the receiver belongs
+  array1d< localIndex > m_receiverRegion;
+
   /// Array containing the elements which contain a receiver
   array1d< localIndex > m_rcvElem;
 

src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp

@@ -66,6 +66,7 @@ struct PrecomputeSourceAndReceiverKernel
   template< typename EXEC_POLICY, typename FE_TYPE >
   static void
   launch( localIndex const size,
+          localIndex const regionIndex,
           localIndex const numNodesPerElem,
           localIndex const numFacesPerElem,
           arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X,
@@ -80,11 +81,13 @@ struct PrecomputeSourceAndReceiverKernel
           arrayView1d< localIndex > const sourceElem,
           arrayView2d< localIndex > const sourceNodeIds,
           arrayView2d< real64 > const sourceConstants,
+          arrayView1d< localIndex > const sourceRegion,
           arrayView2d< real64 const > const receiverCoordinates,
           arrayView1d< localIndex > const receiverIsLocal,
           arrayView1d< localIndex > const rcvElem,
           arrayView2d< localIndex > const receiverNodeIds,
           arrayView2d< real64 > const receiverConstants,
+          arrayView1d< localIndex > const receiverRegion,
           arrayView2d< real32 > const sourceValue,
           real64 const dt,
           real32 const timeSourceFrequency,
@@ -127,6 +130,7 @@ struct PrecomputeSourceAndReceiverKernel
 
             sourceIsAccessible[isrc] = 1;
             sourceElem[isrc] = k;
+            sourceRegion[isrc] = regionIndex;
             real64 Ntest[FE_TYPE::numNodes];
             FE_TYPE::calcN( coordsOnRefElem, Ntest );
 
@@ -174,6 +178,7 @@ struct PrecomputeSourceAndReceiverKernel
                                                                               coordsOnRefElem );
             receiverIsLocal[ircv] = 1;
             rcvElem[ircv] = k;
+            receiverRegion[ircv] = regionIndex;
 
             real64 Ntest[FE_TYPE::numNodes];
             FE_TYPE::calcN( coordsOnRefElem, Ntest );
@@ -443,6 +448,7 @@ struct PressureComputation
    * @tparam EXEC_POLICY the execution policy
    * @tparam ATOMIC_POLICY the atomic policy
    * @param[in] size the number of cells in the subRegion
+   * @param[in] regionIndex Index of the subregion
    * @param[in] size_node the number of nodes in the subRegion
    * @param[in] X coordinates of the nodes
    * @param[in] elemsToNodes map from element to nodes
@@ -463,6 +469,7 @@ struct PressureComputation
   template< typename EXEC_POLICY, typename ATOMIC_POLICY >
   void
   launch( localIndex const size,
+          localIndex const regionIndex,
           localIndex const size_node,
           arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X,
           arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes,
@@ -475,6 +482,7 @@ struct PressureComputation
           arrayView2d< real32 const > const sourceValue,
           arrayView1d< localIndex const > const sourceIsAccessible,
           arrayView1d< localIndex const > const sourceElem,
+          arrayView1d< localIndex const > const sourceRegion,
           real64 const dt,
           integer const cycleNumber,
           arrayView1d< real32 > const p_np1 )
@@ -552,7 +560,7 @@ struct PressureComputation
       {
         if( sourceIsAccessible[isrc] == 1 )
         {
-          if( sourceElem[isrc]==k )
+          if( sourceElem[isrc]==k && sourceRegion[isrc] == regionIndex )
           {
             for( localIndex i = 0; i < numNodesPerElem; ++i )
             {