feat: Add attenuation for acoustic wave solvers

src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFields.hpp

@@ -74,6 +74,22 @@ DECLARE_FIELD( StiffnessVector_q,
                WRITE_AND_READ,
                "Stiffness vector contains R_h*Pressure_n." );
 
+DECLARE_FIELD( StiffnessVectorA_p,
+               "stiffnessVectorA_p",
+               array1d< real32 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "P-type acoustic attenuation stiffness vector." );
+
+DECLARE_FIELD( StiffnessVectorA_q,
+               "stiffnessVectorA_q",
+               array1d< real32 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "Q-type acoustic attenuation stiffness vector." );
+
 DECLARE_FIELD( Pressure_p_nm1,
                "pressure_p_nm1",
                array1d< real32 >,
@@ -122,6 +138,22 @@ DECLARE_FIELD( Pressure_q_np1,
                WRITE_AND_READ,
                "Scalar auxiliary pressure q at time n+1." );
 
+DECLARE_FIELD( DivPsi_p,
+               "divpsi_p",
+               array2d< real32 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "p-type memory variable for acoustic VTI attenuation." );
+
+DECLARE_FIELD( DivPsi_q,
+               "divpsi_q",
+               array2d< real32 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "q-type memory variable for acoustic VTI attenuation." );
+
 DECLARE_FIELD( DampingVector_p,
                "dampingVector_p",
                array1d< real32 >,

src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp

@@ -84,6 +84,16 @@ void AcousticVTIWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
                                acousticvtifields::LateralSurfaceNodeIndicator,
                                acousticvtifields::BottomSurfaceNodeIndicator >( getName() );
 
+    if( m_attenuationType == WaveSolverUtils::AttenuationType::sls )
+    {
+      integer l = m_slsReferenceAngularFrequencies.size( 0 );
+      nodeManager.registerField< acousticvtifields::DivPsi_p,
+                                 acousticvtifields::DivPsi_q,
+                                 acousticvtifields::StiffnessVectorA_p,
+                                 acousticvtifields::StiffnessVectorA_q >( getName() );
+      nodeManager.getField< acousticvtifields::DivPsi_p >().resizeDimension< 1 >( l );
+      nodeManager.getField< acousticvtifields::DivPsi_q >().resizeDimension< 1 >( l );
+    }
 
     FaceManager & faceManager = mesh.getFaceManager();
     faceManager.registerField< acousticfields::AcousticFreeSurfaceFaceIndicator >( getName() );
@@ -98,6 +108,10 @@ void AcousticVTIWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
       subRegion.registerField< acousticvtifields::Epsilon >( getName() );
       subRegion.registerField< acousticvtifields::F >( getName() );
       subRegion.registerField< acousticfields::AcousticVelocity >( getName() );
+      if( m_attenuationType == WaveSolverUtils::AttenuationType::sls )
+      {
+        subRegion.registerField< acousticfields::AcousticQualityFactor >( getName() );
+      }
     } );
   } );
 }
@@ -313,6 +327,12 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
     } );
   } );
 
+  // check anelasticity coefficient and/or compute it if needed
+  if( m_attenuationType == WaveSolverUtils::AttenuationType::sls )
+  {
+    initializeAnelasticityCoefficients< acousticfields::AcousticQualityFactor >( domain );
+  }
+
   WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal );
 }
 
@@ -590,22 +610,56 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n,
                                                             "",
                                                             kernelFactory );
 
+
+    if( m_attenuationType == WaveSolverUtils::AttenuationType::sls )
+    {
+      auto kernelFactory = acousticVTIWaveEquationSEMKernels::ExplicitAcousticVTIAttenuativeSEMFactory( dt );
+      finiteElement::
+        regionBasedKernelApplication< EXEC_POLICY,
+                                      constitutive::NullModel,
+                                      CellElementSubRegion >( mesh,
+                                                              regionNames,
+                                                              getDiscretizationName(),
+                                                              "",
+                                                              kernelFactory );
+    }
+
     addSourceToRightHandSide( time_n, rhs );
 
     /// calculate your time integrators
 
     GEOS_MARK_SCOPE ( updateP );
 
-    AcousticTimeSchemeSEM::LeapFrogforVTI( nodeManager.size(), dt, p_np1, p_n, p_nm1, q_np1, q_n, q_nm1, mass, stiffnessVector_p,
-                                           stiffnessVector_q, damping_p, damping_pq, damping_q, damping_qp,
-                                           rhs, freeSurfaceNodeIndicator, lateralSurfaceNodeIndicator,
-                                           bottomSurfaceNodeIndicator );
-
+    if( m_attenuationType == WaveSolverUtils::AttenuationType::sls )
+    {
+      arrayView1d< real32 > const stiffnessVectorA_p = nodeManager.getField< acousticvtifields::StiffnessVectorA_p >();
+      arrayView1d< real32 > const stiffnessVectorA_q = nodeManager.getField< acousticvtifields::StiffnessVectorA_q >();
+      arrayView2d< real32 > const divpsi_p = nodeManager.getField< acousticvtifields::DivPsi_p >();
+      arrayView2d< real32 > const divpsi_q = nodeManager.getField< acousticvtifields::DivPsi_q >();
+      arrayView1d< real32 > const referenceFrequencies = m_slsReferenceAngularFrequencies.toView();
+      arrayView1d< real32 > const anelasticityCoefficients = m_slsAnelasticityCoefficients.toView();
+      AcousticTimeSchemeSEM::AttenuationLeapFrogforVTI( nodeManager.size(), dt, p_np1, p_n, p_nm1, q_np1, q_n, q_nm1, divpsi_p, divpsi_q, mass, stiffnessVector_p,
+                                                        stiffnessVector_q, stiffnessVectorA_p, stiffnessVectorA_q, damping_p, damping_pq, damping_q, damping_qp, rhs, 
+                                                        freeSurfaceNodeIndicator, lateralSurfaceNodeIndicator, bottomSurfaceNodeIndicator, referenceFrequencies,
+                                                        anelasticityCoefficients );
+    }
+    else
+    {
+      AcousticTimeSchemeSEM::LeapFrogforVTI( nodeManager.size(), dt, p_np1, p_n, p_nm1, q_np1, q_n, q_nm1, mass, stiffnessVector_p,
+                                             stiffnessVector_q, damping_p, damping_pq, damping_q, damping_qp,
+                                             rhs, freeSurfaceNodeIndicator, lateralSurfaceNodeIndicator,
+                                             bottomSurfaceNodeIndicator );
+    }
     /// synchronize pressure fields
     FieldIdentifiers fieldsToBeSync;
     fieldsToBeSync.addFields( FieldLocation::Node, { acousticvtifields::Pressure_p_np1::key() } );
     fieldsToBeSync.addFields( FieldLocation::Node, { acousticvtifields::Pressure_q_np1::key() } );
 
+   if( m_slsReferenceAngularFrequencies.size( 0 ) > 0 )
+   {
+     fieldsToBeSync.addFields( FieldLocation::Node, { acousticvtifields::DivPsi_p::key(), acousticvtifields::DivPsi_q::key() } );
+   }
+
     CommunicationTools & syncFields = CommunicationTools::getInstance();
     syncFields.synchronizeFields( fieldsToBeSync,
                                   mesh,
@@ -625,6 +679,15 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n,
       stiffnessVector_q[a] = 0.0;
       rhs[a] = 0.0;
     } );
+    if( m_attenuationType == WaveSolverUtils::AttenuationType::sls )
+    {
+      arrayView1d< real32 > const stiffnessVectorA_p = nodeManager.getField< acousticvtifields::StiffnessVectorA_p >();
+      arrayView1d< real32 > const stiffnessVectorA_q = nodeManager.getField< acousticvtifields::StiffnessVectorA_q >();
+      forAll< EXEC_POLICY >( nodeManager.size() , [=] GEOS_HOST_DEVICE ( localIndex const a )
+      {
+        stiffnessVectorA_p[a] = stiffnessVectorA_q[a] = 0.0;
+      } );
+    }
 
   } );
 

src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp

@@ -376,12 +376,13 @@ struct DampingMatrixKernel
 
 template< typename SUBREGION_TYPE,
           typename CONSTITUTIVE_TYPE,
-          typename FE_TYPE >
-class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
-                                                                 CONSTITUTIVE_TYPE,
-                                                                 FE_TYPE,
-                                                                 1,
-                                                                 1 >
+          typename FE_TYPE,
+          typename Sp = fields::acousticvtifields::StiffnessVector_p, typename Sq = fields::acousticvtifields::StiffnessVector_q >
+class ExplicitAcousticVTISEMBase : public finiteElement::KernelBase< SUBREGION_TYPE,
+                                                                     CONSTITUTIVE_TYPE,
+                                                                     FE_TYPE,
+                                                                     1,
+                                                                     1 >
 {
 public:
 
@@ -415,22 +416,22 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
    * @param dt The time interval for the step.
    *   elements to be processed during this kernel launch.
    */
-  ExplicitAcousticVTISEM( NodeManager & nodeManager,
-                          EdgeManager const & edgeManager,
-                          FaceManager const & faceManager,
-                          localIndex const targetRegionIndex,
-                          SUBREGION_TYPE const & elementSubRegion,
-                          FE_TYPE const & finiteElementSpace,
-                          CONSTITUTIVE_TYPE & inputConstitutiveType,
-                          real64 const dt ):
+  ExplicitAcousticVTISEMBase( NodeManager & nodeManager,
+                              EdgeManager const & edgeManager,
+                              FaceManager const & faceManager,
+                              localIndex const targetRegionIndex,
+                              SUBREGION_TYPE const & elementSubRegion,
+                              FE_TYPE const & finiteElementSpace,
+                              CONSTITUTIVE_TYPE & inputConstitutiveType,
+                              real64 const dt ):
     Base( elementSubRegion,
           finiteElementSpace,
           inputConstitutiveType ),
     m_nodeCoords( nodeManager.getField< fields::referencePosition32 >() ),
     m_p_n( nodeManager.getField< fields::acousticvtifields::Pressure_p_n >() ),
     m_q_n( nodeManager.getField< fields::acousticvtifields::Pressure_q_n >() ),
-    m_stiffnessVector_p( nodeManager.getField< fields::acousticvtifields::StiffnessVector_p >() ),
-    m_stiffnessVector_q( nodeManager.getField< fields::acousticvtifields::StiffnessVector_q >() ),
+    m_stiffnessVector_p( nodeManager.getField< Sp >() ),
+    m_stiffnessVector_q( nodeManager.getField< Sq >() ),
     m_epsilon( elementSubRegion.template getField< fields::acousticvtifields::Epsilon >() ),
     m_delta( elementSubRegion.template getField< fields::acousticvtifields::Delta >() ),
     m_vti_f( elementSubRegion.template getField< fields::acousticvtifields::F >() ),
@@ -445,7 +446,7 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
   /**
    * @copydoc geos::finiteElement::KernelBase::StackVariables
    *
-   * ### ExplicitAcousticVTISEM Description
+   * ### ExplicitAcousticVTISEMBase Description
    * Adds a stack arrays for the nodal force, primary displacement variable, etc.
    */
   struct StackVariables : Base::StackVariables
@@ -461,6 +462,7 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
     /// C-array stack storage for element local the nodal positions.
     /// only the eight corners of the mesh cell are needed to compute the Jacobian
     real64 xLocal[ 8 ][ 3 ];
+    real32 factor;
     /// local (to this element) stiffness vectors
     real32 stiffnessVectorLocal_p[ numNodesPerElem ]{};
     real32 stiffnessVectorLocal_q[ numNodesPerElem ]{};
@@ -487,6 +489,7 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
         stack.xLocal[ a ][ i ] = m_nodeCoords[ nodeIndex ][ i ];
       }
     }
+    stack.factor = 1.0;
   }
 
   GEOS_HOST_DEVICE
@@ -504,7 +507,7 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
   /**
    * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel
    *
-   * ### ExplicitAcousticVTISEM Description
+   * ### ExplicitAcousticVTISEMBase Description
    * Calculates stiffness vector
    *
    */
@@ -518,20 +521,21 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
     m_finiteElementSpace.template computeStiffnessxyTerm( q, stack.xLocal, [&] ( int i, int j, real64 val )
     {
       real32 const localIncrement_p = val*(-1-2*m_epsilon[k])*m_p_n[m_elemsToNodes[k][j]];
-      stack.stiffnessVectorLocal_p[ i ] += localIncrement_p;
+      stack.stiffnessVectorLocal_p[ i ] += localIncrement_p * stack.factor;
+
       real32 const localIncrement_q = val*((-2*m_delta[k]-m_vti_f[k])*m_p_n[m_elemsToNodes[k][j]] +(m_vti_f[k]-1)*m_q_n[m_elemsToNodes[k][j]]);
-      stack.stiffnessVectorLocal_q[ i ] += localIncrement_q;
+      stack.stiffnessVectorLocal_q[ i ] += localIncrement_q * stack.factor;
     } );
 
     // Pseudo-Stiffness z
 
     m_finiteElementSpace.template computeStiffnesszTerm( q, stack.xLocal, [&] ( int i, int j, real64 val )
     {
       real32 const localIncrement_p = val*((m_vti_f[k]-1)*m_p_n[m_elemsToNodes[k][j]] - m_vti_f[k]*m_q_n[m_elemsToNodes[k][j]]);
-      stack.stiffnessVectorLocal_p[ i ] += localIncrement_p;
+      stack.stiffnessVectorLocal_p[ i ] += localIncrement_p * stack.factor;
 
       real32 const localIncrement_q = -val*m_q_n[m_elemsToNodes[k][j]];
-      stack.stiffnessVectorLocal_q[ i ] += localIncrement_q;
+      stack.stiffnessVectorLocal_q[ i ] += localIncrement_q * stack.factor;
     } );
   }
 
@@ -567,10 +571,84 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE,
 };
 
 
-
-/// The factory used to construct a ExplicitAcousticWaveEquation kernel.
+/// Specialization for standard iso elastic kernel
+template< typename SUBREGION_TYPE,
+          typename CONSTITUTIVE_TYPE,
+          typename FE_TYPE >
+using ExplicitAcousticVTISEM = ExplicitAcousticVTISEMBase< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >;
 using ExplicitAcousticVTISEMFactory = finiteElement::KernelFactory< ExplicitAcousticVTISEM,
                                                                     real64 >;
+/// Specialization for attenuation kernel
+template< typename SUBREGION_TYPE,
+          typename CONSTITUTIVE_TYPE,
+          typename FE_TYPE >
+class ExplicitAcousticVTIAttenuativeSEM : public ExplicitAcousticVTISEMBase< SUBREGION_TYPE,
+                                                                             CONSTITUTIVE_TYPE,
+                                                                             FE_TYPE,
+                                                                             fields::acousticvtifields::StiffnessVectorA_p,
+                                                                             fields::acousticvtifields::StiffnessVectorA_q >
+{
+public:
+
+  /// Alias for the base class;
+  using Base = ExplicitAcousticVTISEMBase< SUBREGION_TYPE,
+                                           CONSTITUTIVE_TYPE,
+                                           FE_TYPE,
+                                           fields::acousticvtifields::StiffnessVectorA_p,
+                                           fields::acousticvtifields::StiffnessVectorA_q >;
+
+//*****************************************************************************
+  /**
+   * @brief Constructor
+   * @copydoc geos::finiteElement::KernelBase::KernelBase
+   * @param nodeManager Reference to the NodeManager object.
+   * @param edgeManager Reference to the EdgeManager object.
+   * @param faceManager Reference to the FaceManager object.
+   * @param targetRegionIndex Index of the region the subregion belongs to.
+   * @param dt The time interval for the step.
+   */
+  ExplicitAcousticVTIAttenuativeSEM( NodeManager & nodeManager,
+                                     EdgeManager const & edgeManager,
+                                     FaceManager const & faceManager,
+                                     localIndex const targetRegionIndex,
+                                     SUBREGION_TYPE const & elementSubRegion,
+                                     FE_TYPE const & finiteElementSpace,
+                                     CONSTITUTIVE_TYPE & inputConstitutiveType,
+                                     real64 const dt ):
+    Base( nodeManager,
+          edgeManager,
+          faceManager,
+          targetRegionIndex,
+          elementSubRegion,
+          finiteElementSpace,
+          inputConstitutiveType,
+          dt ),
+    m_qualityFactor( elementSubRegion.template getField< fields::acousticfields::AcousticQualityFactor >() )
+  {}
+
+  /**
+   * @copydoc geos::finiteElement::KernelBase::setup
+   *
+   * Copies the primary variable, and position into the local stack array.
+   */
+  GEOS_HOST_DEVICE
+  inline
+  void setup( localIndex const k,
+              typename Base::StackVariables & stack ) const
+  {
+    Base::setup( k, stack );
+    stack.factor = stack.factor / m_qualityFactor[ k ];
+  }
+
+protected:
+
+  /// The array containing the acoustic attenuation quality factor
+  arrayView1d< real32 const > const m_qualityFactor;
+
+};
+
+using ExplicitAcousticVTIAttenuativeSEMFactory = finiteElement::KernelFactory< ExplicitAcousticVTIAttenuativeSEM,
+                                                                               real64 >;
 
 } // namespace acousticVTIWaveEquationSEMKernels