rewrite seismo traces - rework WaveSolverBase structure

src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp

@@ -65,21 +65,10 @@ AcousticFirstOrderWaveEquationSEM::AcousticFirstOrderWaveEquationSEM( const std:
     setSizedFromParent( 0 ).
     setDescription( "Element containing the sources" );
 
-  registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ).
-    setInputFlag( InputFlags::FALSE ).
-    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()
@@ -151,10 +140,10 @@ void AcousticFirstOrderWaveEquationSEM::postProcessInput()
 
   localIndex const numSourcesGlobal = m_sourceCoordinates.size( 0 );
   localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 );
-  m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
-  m_uxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
-  m_uyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
-  m_uzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
+  m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+  m_uxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+  m_uyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+  m_uzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
   m_sourceElem.resize( numSourcesGlobal );
   m_sourceRegion.resize( numSourcesGlobal );
   m_rcvElem.resize( numReceiversGlobal );
@@ -288,8 +277,7 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
 
   DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
 
-  real64 const time = 0.0;
-  applyFreeSurfaceBC( time, domain );
+  applyFreeSurfaceBC( 0.0, domain );
 
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                 MeshLevel & mesh,
@@ -354,6 +342,8 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
       } );
     } );
   } );
+
+  WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
 }
 
 
@@ -524,9 +514,9 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
           cycleNumber,
           p_np1 );
       } );
-      arrayView2d< real32 > const uxReceivers   = m_uxNp1AtReceivers.toView();
-      arrayView2d< real32 > const uyReceivers   = m_uyNp1AtReceivers.toView();
-      arrayView2d< real32 > const uzReceivers   = m_uzNp1AtReceivers.toView();
+      arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView();
+      arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView();
+      arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView();
 
       compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_x, velocity_x, uxReceivers );
       compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_y, velocity_y, uyReceivers );
@@ -545,16 +535,12 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
                                   true );
 
     // compute the seismic traces since last step.
-    arrayView2d< real32 > const pReceivers   = m_pressureNp1AtReceivers.toView();
+    arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
     computeAllSeismoTraces( time_n, dt, p_np1, p_np1, pReceivers );
 
-    // increment m_indexSeismoTrace
-    while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
-    {
-      m_indexSeismoTrace++;
-    }
-
+    incrementIndexSeismoTrace( time_n );
   } );
+
   return dt;
 }
 
@@ -574,62 +560,26 @@ void AcousticFirstOrderWaveEquationSEM::cleanup( real64 const time_n, integer co
       arrayView2d< real32 > const velocity_y = elementSubRegion.getField< wavesolverfields::Velocity_y >();
       arrayView2d< real32 > const velocity_z = elementSubRegion.getField< wavesolverfields::Velocity_z >();
 
-      arrayView2d< real32 > const uxReceivers   = m_uxNp1AtReceivers.toView();
-      arrayView2d< real32 > const uyReceivers   = m_uyNp1AtReceivers.toView();
-      arrayView2d< real32 > const uzReceivers   = m_uzNp1AtReceivers.toView();
-
-      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 uxReceivers = m_uxNp1AtReceivers.toView();
+      arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView();
+      arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView();
 
-    } );
-    arrayView2d< real32 > const pReceivers   = m_pressureNp1AtReceivers.toView();
-    computeAllSeismoTraces( time_n, 0, p_np1, p_np1, pReceivers );
-  } );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_x, velocity_x, uxReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_y, velocity_y, uyReceivers );
+      compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_z, velocity_z, uzReceivers );
 
-  // increment m_indexSeismoTrace
-  while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
-  {
-    m_indexSeismoTrace++;
-  }
-}
+      WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+                                               m_receiverIsLocal, m_nsamplesSeismoTrace, uxReceivers, uyReceivers, uzReceivers );
 
-void AcousticFirstOrderWaveEquationSEM::computeAllSeismoTraces( real64 const time_n,
-                                                                real64 const dt,
-                                                                arrayView1d< real32 const > const var_np1,
-                                                                arrayView1d< real32 const > const var_n,
-                                                                arrayView2d< real32 > varAtReceivers )
-{
-  localIndex indexSeismoTrace = m_indexSeismoTrace;
-  for( real64 timeSeismo;
-       (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
-       indexSeismoTrace++ )
-  {
-    WaveSolverUtils::computeSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal, m_nsamplesSeismoTrace,
-                                         m_outputSeismoTrace,
-                                         var_np1, var_n, varAtReceivers );
-  }
-}
+    } );
+    arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
+    computeAllSeismoTraces( time_n, 0.0, p_np1, p_np1, pReceivers );
+    WaveSolverUtils::writeSeismoTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+                                       m_receiverIsLocal, m_nsamplesSeismoTrace, pReceivers );
 
-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,
-                                                                          arrayView2d< real32 > varAtReceivers )
-{
-  localIndex indexSeismoTrace = m_indexSeismoTrace;
-  for( real64 timeSeismo;
-       (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
-       indexSeismoTrace++ )
-  {
-    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 );
-  }
+  } );
 }
 
-
 void AcousticFirstOrderWaveEquationSEM::initializePML()
 {
   GEOS_ERROR( getDataContext() << ": PML for the first order acoustic wave propagator not yet implemented" );

src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp

@@ -21,8 +21,8 @@
 #define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICFIRSTORDERWAVEEQUATIONSEM_HPP_
 
 #include "mesh/MeshFields.hpp"
-#include "WaveSolverUtils.hpp"
 #include "WaveSolverBaseFields.hpp"
+#include "WaveSolverBase.hpp"
 
 namespace geos
 {
@@ -78,37 +78,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
    */
   virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs );
 
-  /**
-   * TODO: move implementation into WaveSolverUtils
-   * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
-   * @param time_n the time corresponding to the field values pressure_n
-   * @param dt the simulation timestep
-   * @param var_np1 the field values at time_n + dt
-   * @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 computeAllSeismoTraces( real64 const time_n,
-                                       real64 const dt,
-                                       arrayView1d< real32 const > const var_np1,
-                                       arrayView1d< real32 const > const var_n,
-                                       arrayView2d< real32 > varAtReceivers );
-
-  /**
-   * TODO: move implementation into WaveSolverUtils
-   * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt for a 2d variable
-   * @param time_n the time corresponding to the field values pressure_n
-   * @param dt the simulation timestep
-   * @param var_np1 the field values at time_n + dt
-   * @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( localIndex const regionIndex,
-                                                 real64 const time_n,
-                                                 real64 const dt,
-                                                 arrayView2d< real32 const > const var_np1,
-                                                 arrayView2d< real32 const > const var_n,
-                                                 arrayView2d< real32 > varAtReceivers );
-
 
   /**
    * @brief Initialize Perfectly Matched Layer (PML) information
@@ -133,8 +102,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
 
     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;
 
@@ -197,13 +164,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
 
   /// 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/AcousticVTIWaveEquationSEM.cpp

@@ -106,7 +106,7 @@ void AcousticVTIWaveEquationSEM::postProcessInput()
 
   localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 );
 
-  m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
+  m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
 }
 
 void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
@@ -229,8 +229,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
 
   DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
 
-  real64 const time = 0.0;
-  applyFreeSurfaceBC( time, domain );
+  applyFreeSurfaceBC( 0.0, domain );
   precomputeSurfaceFieldIndicator( domain );
 
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
@@ -315,6 +314,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
     } );
   } );
 
+  WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
 }
 
 void AcousticVTIWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain )
@@ -649,10 +649,11 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n,
                                   true );
 
     // compute the seismic traces since last step.
-    arrayView2d< real32 > const pReceivers   = m_pressureNp1AtReceivers.toView();
-
+    arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
     computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers );
 
+    incrementIndexSeismoTrace( time_n );
+
     /// prepare next step
     forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
     {
@@ -681,45 +682,16 @@ void AcousticVTIWaveEquationSEM::cleanup( real64 const time_n,
                                                                 arrayView1d< string const > const & )
   {
     NodeManager & nodeManager = mesh.getNodeManager();
-    arrayView1d< real32 const > const p_n = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >();
+    arrayView1d< real32 const > const p_n   = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >();
     arrayView1d< real32 const > const p_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_np1 >();
-    arrayView2d< real32 > const pReceivers   = m_pressureNp1AtReceivers.toView();
-    computeAllSeismoTraces( time_n, 0, p_np1, p_n, pReceivers );
-  } );
-}
+    arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
+    computeAllSeismoTraces( time_n, 0.0, p_np1, p_n, pReceivers );
 
-void AcousticVTIWaveEquationSEM::computeAllSeismoTraces( real64 const time_n,
-                                                         real64 const dt,
-                                                         arrayView1d< real32 const > const var_np1,
-                                                         arrayView1d< real32 const > const var_n,
-                                                         arrayView2d< real32 > varAtReceivers )
-{
-
-  /*
-   * In forward case we compute seismo if time_n + dt  is the first time
-   * step after the timeSeismo to write.
-   *
-   *  time_n        timeSeismo    time_n + dt
-   *   ---|--------------|-------------|
-   *
-   * In backward (time_n goes decreasing) case we compute seismo if
-   * time_n is the last time step before the timeSeismo to write.
-   *
-   *  time_n - dt    timeSeismo    time_n
-   *   ---|--------------|-------------|
-   */
-  for( real64 timeSeismo;
-       (m_forward)?((timeSeismo = m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + dt + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace):
-       ((timeSeismo = m_dtSeismoTrace*(m_nsamplesSeismoTrace-m_indexSeismoTrace-1)) >= (time_n - dt -  epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace);
-       m_indexSeismoTrace++ )
-  {
-    WaveSolverUtils::computeSeismoTrace( time_n, (m_forward)?dt:-dt, timeSeismo, (m_forward)?m_indexSeismoTrace:(m_nsamplesSeismoTrace-m_indexSeismoTrace-1), m_receiverNodeIds, m_receiverConstants,
-                                         m_receiverIsLocal,
-                                         m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers );
-  }
+    WaveSolverUtils::writeSeismoTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+                                       m_receiverIsLocal, m_nsamplesSeismoTrace, pReceivers );
+  } );
 }
 
-
 REGISTER_CATALOG_ENTRY( SolverBase, AcousticVTIWaveEquationSEM, string const &, dataRepository::Group * const )
 
 } /* namespace geos */

src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp

@@ -72,21 +72,6 @@ class AcousticVTIWaveEquationSEM : public WaveSolverBase
    */
   virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs );
 
-  /**
-   * TODO: move implementation into WaveSolverBase
-   * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
-   * @param time_n the time corresponding to the field values pressure_n
-   * @param dt the simulation timestep
-   * @param var_np1 the field values at time_n + dt
-   * @param var_n the field values at time_n
-   * @param varAtReceivers the array holding the trace values, where the output is written
-   */
-  virtual void computeAllSeismoTraces( real64 const time_n,
-                                       real64 const dt,
-                                       arrayView1d< real32 const > const var_np1,
-                                       arrayView1d< real32 const > const var_n,
-                                       arrayView2d< real32 > varAtReceivers );
-
   /**
    * @brief Overridden from ExecutableGroup. Used to write last seismogram if needed.
    */