diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFields.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFields.hpp index a53cf0c52a..5fd3c040aa 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFields.hpp +++ b/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 >, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp index 28b2497012..c0afb1a484 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +++ b/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; + } ); + } } ); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp index c0b6e66e19..88de22a519 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp +++ b/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,9 +521,10 @@ 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 @@ -528,10 +532,10 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE, 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 diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp index c47e6e83fc..c1dd685e9f 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -84,6 +84,14 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) acousticfields::StiffnessVector, acousticfields::AcousticFreeSurfaceNodeIndicator >( getName() ); + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + integer l = m_slsReferenceAngularFrequencies.size( 0 ); + nodeManager.registerField< acousticfields::DivPsi, + acousticfields::StiffnessVectorA >( getName() ); + nodeManager.getField< acousticfields::DivPsi >().resizeDimension< 1 >( l ); + } + /// register PML auxiliary variables only when a PML is specified in the xml if( m_usePML ) { @@ -106,6 +114,10 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< acousticfields::AcousticVelocity >( getName() ); subRegion.registerField< acousticfields::AcousticDensity >( getName() ); subRegion.registerField< acousticfields::PartialGradient >( getName() ); + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + subRegion.registerField< acousticfields::AcousticQualityFactor >( getName() ); + } subRegion.registerField< acousticfields::PartialGradient2 >( getName() ); } ); @@ -324,12 +336,16 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() velocity, density, damping ); - - } ); } ); } ); + // check anelasticity coefficient and/or compute it if needed + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + initializeAnelasticityCoefficients< acousticfields::AcousticQualityFactor >( domain ); + } + if( m_timestepStabilityLimit==1 ) { real64 dtOut = 0.0; @@ -1072,6 +1088,15 @@ void AcousticWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh ) stiffnessVector[a] = rhs[a] = 0.0; } ); + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< acousticfields::StiffnessVectorA >(); + forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n ) + { + localIndex const a = solverTargetNodesSet[n]; + stiffnessVectorA[a] = 0.0; + } ); + } } void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, @@ -1103,6 +1128,18 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, getDiscretizationName(), "", kernelFactory ); + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + auto kernelFactory = acousticWaveEquationSEMKernels::ExplicitAcousticAttenuativeSEMFactory( dt ); + finiteElement:: + regionBasedKernelApplication< EXEC_POLICY, + constitutive::NullModel, + CellElementSubRegion >( mesh, + regionNames, + getDiscretizationName(), + "", + kernelFactory ); + } //Modification of cycleNember useful when minTime < 0 addSourceToRightHandSide( time_n, rhs ); @@ -1110,11 +1147,28 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, real64 const dt2 = pow( dt, 2 ); SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst(); + if( m_usePML && m_attenuationType != WaveSolverUtils::AttenuationType::none ) + { + GEOS_ERROR( "Attenuation is not supported with PML boindary conditions."); + } if( !m_usePML ) { GEOS_MARK_SCOPE ( updateP ); - AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, - rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< acousticfields::StiffnessVectorA >(); + arrayView2d< real32 > const divpsi = nodeManager.getField< acousticfields::DivPsi >(); + arrayView1d< real32 > const referenceFrequencies = m_slsReferenceAngularFrequencies.toView(); + arrayView1d< real32 > const anelasticityCoefficients = m_slsAnelasticityCoefficients.toView(); + AcousticTimeSchemeSEM::AttenuationLeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, divpsi, mass, stiffnessVector, stiffnessVectorA, + damping, rhs, freeSurfaceNodeIndicator, solverTargetNodesSet, referenceFrequencies, + anelasticityCoefficients ); + } + else + { + AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, + rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); + } } else { @@ -1198,6 +1252,11 @@ void AcousticWaveEquationSEM::synchronizeUnknowns( real64 const & time_n, FieldIdentifiers fieldsToBeSync; fieldsToBeSync.addFields( FieldLocation::Node, { acousticfields::Pressure_np1::key() } ); + if( m_slsReferenceAngularFrequencies.size( 0 ) > 0 ) + { + fieldsToBeSync.addFields( FieldLocation::Node, { acousticfields::DivPsi::key() } ); + } + if( m_usePML ) { fieldsToBeSync.addFields( FieldLocation::Node, { diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEMKernel.hpp index 1e993b9aa4..89a9ea4f10 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEMKernel.hpp @@ -53,12 +53,13 @@ namespace acousticWaveEquationSEMKernels template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, - typename FE_TYPE > -class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, - CONSTITUTIVE_TYPE, - FE_TYPE, - 1, - 1 > + typename FE_TYPE, + typename S = fields::acousticfields::StiffnessVector > +class ExplicitAcousticSEMBase : public finiteElement::KernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 1, + 1 > { public: @@ -92,20 +93,20 @@ class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, * @param dt The time interval for the step. * elements to be processed during this kernel launch. */ - ExplicitAcousticSEM( 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 ): + ExplicitAcousticSEMBase( 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::acousticfields::Pressure_n >() ), - m_stiffnessVector( nodeManager.getField< fields::acousticfields::StiffnessVector >() ), + m_stiffnessVector( nodeManager.getField< S >() ), m_density( elementSubRegion.template getField< fields::acousticfields::AcousticDensity >() ), m_dt( dt ) { @@ -118,7 +119,7 @@ class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, /** * @copydoc geos::finiteElement::KernelBase::StackVariables * - * ### ExplicitAcousticSEM Description + * ### ExplicitAcousticSEMBase Description * Adds a stack arrays for the nodal force, primary displacement variable, etc. */ struct StackVariables : Base::StackVariables @@ -178,7 +179,7 @@ class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, /** * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel * - * ### ExplicitAcousticSEM Description + * ### ExplicitAcousticSEMBase Description * Calculates stiffness vector * */ @@ -217,10 +218,84 @@ class ExplicitAcousticSEM : 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 ExplicitAcousticSEM = ExplicitAcousticSEMBase< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >; using ExplicitAcousticSEMFactory = finiteElement::KernelFactory< ExplicitAcousticSEM, real64 >; +/// Specialization for attenuation kernel +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class ExplicitAcousticAttenuativeSEM : public ExplicitAcousticSEMBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + fields::acousticfields::StiffnessVectorA > +{ +public: + + /// Alias for the base class; + using Base = ExplicitAcousticSEMBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + fields::acousticfields::StiffnessVectorA >; + +//***************************************************************************** + /** + * @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. + */ + ExplicitAcousticAttenuativeSEM( 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.invDensity = stack.invDensity / m_qualityFactor[ k ]; + } + +protected: + + /// The array containing the acoustic attenuation quality factor + arrayView1d< real32 const > const m_qualityFactor; + +}; + +using ExplicitAcousticAttenuativeSEMFactory = finiteElement::KernelFactory< ExplicitAcousticAttenuativeSEM, + real64 >; + } // namespace acousticWaveEquationSEMKernels diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp index 5222bc115a..c220d3063a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp @@ -59,6 +59,30 @@ DECLARE_FIELD( Pressure_np1, WRITE_AND_READ, "Scalar pressure at time n+1." ); +DECLARE_FIELD( DivPsi, + "divpsi", + array2d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "memory variable for acoustic attenuation." ); + +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( PressureForward, "pressureForward", array1d< real32 >, @@ -131,6 +155,14 @@ DECLARE_FIELD( StiffnessVector, WRITE_AND_READ, "Stiffness vector contains R_h*Pressure_n." ); +DECLARE_FIELD( StiffnessVectorA, + "stiffnessVectorA", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Acoustic attenuation stiffness vector." ); + DECLARE_FIELD( DampingVector, "dampingVector", array1d< real32 >, @@ -155,6 +187,14 @@ DECLARE_FIELD( AcousticDensity, WRITE_AND_READ, "Medium density of the cell" ); +DECLARE_FIELD( AcousticQualityFactor, + "acousticQualityFactor", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Quality factor for acoustic wave attenuation in the cell" ); + DECLARE_FIELD( AcousticFreeSurfaceFaceIndicator, "acousticFreeSurfaceFaceIndicator", array1d< localIndex >, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp index 69b78be7a7..1c90dc2fe4 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp @@ -69,6 +69,64 @@ struct AcousticTimeSchemeSEM }; + /** + * @brief Apply second order Leap-Frog time scheme for isotropic case without PML, but with attenuation + * @param[in] dt time-step + * @param[out] p_np1 pressure array at time n+1 (updated here) + * @param[in] p_n pressure array at time n + * @param[in] p_nm1 pressure array at time n-1 + * @param[in] divpsi divergence of the acoustic memory variable + * @param[in] mass the mass matrix + * @param[in] stiffnessVector array containing the product of the stiffness matrix R and the pressure at time n + * @param[in] stiffnessVectorA array containing the product of the attenuation stiffness matrix R and the pressure at time n + * @param[in] damping the damping matrix + * @param[in] rhs the right-hand-side + * @param[in] freeSurfaceNodeIndicator array which contains indicators to tell if we are on a free-surface boundary or not + * @param[in] solverTargetNodesSet the targetted nodeset (useful in particular when we do elasto-acoustic simulation ) + * @param[in] referenceFrequencies the reference frequencies for each SLS + * @param[in] anelasticityCoefficients the coefficients of anelasticity for each SLS + */ + static void AttenuationLeapFrogWithoutPML( real64 const dt, + arrayView1d< real32 > const p_np1, + arrayView1d< real32 > const p_n, + arrayView1d< real32 > const p_nm1, + arrayView2d< real32 > const divpsi, + arrayView1d< real32 const > const mass, + arrayView1d< real32 > const stiffnessVector, + arrayView1d< real32 > const stiffnessVectorA, + arrayView1d< real32 const > const damping, + arrayView1d< real32 > const rhs, + arrayView1d< localIndex const > const freeSurfaceNodeIndicator, + SortedArrayView< localIndex const > const solverTargetNodesSet, + arrayView1d< real32 > referenceFrequencies, + arrayView1d< real32 > anelasticityCoefficients ) + { + real64 const dt2 = pow( dt, 2 ); + integer nl = referenceFrequencies.size( 0 ); + forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n ) + { + localIndex const a = solverTargetNodesSet[n]; + if( freeSurfaceNodeIndicator[a] != 1 ) + { + p_np1[a] = p_n[a]; + p_np1[a] *= 2.0 * mass[a]; + p_np1[a] -= (mass[a] - 0.5 * dt * damping[a]) * p_nm1[a]; + p_np1[a] += dt2 * (rhs[a] - stiffnessVector[a]); + for( integer l = 0; l < nl; l++ ) + { + real32 gammal = ( 2.0 - referenceFrequencies[ l ] * dt )/( 2.0 + referenceFrequencies[ l ] * dt ); + real32 betal = anelasticityCoefficients[ l ] * referenceFrequencies[ l ] * 2.0 * dt /( 2.0 + referenceFrequencies[ l ] * dt ); + real32 gammalp = 0.5 + 0.5 * gammal; + real32 betalp = 0.5 * betal; + p_np1[a] += dt2 * ( gammalp * divpsi( a, l ) + betalp * stiffnessVectorA[ a ] ); + divpsi( a, l ) = gammal * divpsi( a, l ) + betal * stiffnessVectorA[ a ]; + } + p_np1[a] /= mass[a] + 0.5 * dt * damping[a]; + } + } ); + + }; + /** * @brief Apply second order Leap-Frog time scheme for VTI case without PML * @param[in] size The number of nodes in the nodeManager @@ -160,6 +218,122 @@ struct AcousticTimeSchemeSEM } ); }; + /** + * @brief Apply second order Leap-Frog time scheme for VTI case without PML, with attenuation + * @param[in] size The number of nodes in the nodeManager + * @param[in] dt time-step + * @param[out] p_np1 pressure array at time n+1 (updated here) + * @param[in] p_n pressure array at time n + * @param[in] p_nm1 pressure array at time n-1 + * @param[out] q_np1 auxiliary pressure array at time n+1 (updated here) + * @param[in] q_n auxiliary pressure array at time n + * @param[in] q_nm1 auxiliary pressure array at time n-1 + * @param[in] divpsi_p divergence of the p-type acoustc memory variable + * @param[in] divpsi_q divergence of the q-type acoustc memory variable + * @param[in] mass the mass matrix + * @param[in] stiffnessVector_p array containing the product of the stiffness matrix R and the pressure at time n + * @param[in] stiffnessVector_q array containing the product of the stiffness matrix R and the auxiliary pressure at time n + * @param[in] stiffnessVectorA_p array containing the product of the attenuated stiffness matrix R and the pressure at time n + * @param[in] stiffnessVectorA_q array containing the product of the attenuated stiffness matrix R and the auxiliary pressure at time n + * @param[in] damping_p the damping matrix + * @param[in] damping_pq the damping matrix + * @param[in] damping_q the damping matrix + * @param[in] damping_qp the damping matrix + * @param[in] rhs the right-hand-side + * @param[in] freeSurfaceNodeIndicator array which contains indicators to tell if we are on a free-surface boundary or not + * @param[in] lateralSurfaceNodeIndicator array which contains indicators to tell if we are on a lateral boundary or not + * @param[in] bottomSurfaceNodeIndicator array which contains indicators to telle if we are on the bottom boundary or not + * @param[in] referenceFrequencies the reference frequencies for each SLS + * @param[in] anelasticityCoefficients the coefficients of anelasticity for each SLS + */ + static void AttenuationLeapFrogforVTI( localIndex const size, + real64 const dt, + arrayView1d< real32 > const p_np1, + arrayView1d< real32 > const p_n, + arrayView1d< real32 > const p_nm1, + arrayView1d< real32 > const q_np1, + arrayView1d< real32 > const q_n, + arrayView1d< real32 > const q_nm1, + arrayView2d< real32 > const divpsi_p, + arrayView2d< real32 > const divpsi_q, + arrayView1d< real32 const > const mass, + arrayView1d< real32 > const stiffnessVector_p, + arrayView1d< real32 > const stiffnessVector_q, + arrayView1d< real32 > const stiffnessVectorA_p, + arrayView1d< real32 > const stiffnessVectorA_q, + arrayView1d< real32 const > const damping_p, + arrayView1d< real32 const > const damping_pq, + arrayView1d< real32 const > const damping_q, + arrayView1d< real32 const > const damping_qp, + arrayView1d< real32 > const rhs, + arrayView1d< localIndex const > const freeSurfaceNodeIndicator, + arrayView1d< localIndex const > const lateralSurfaceNodeIndicator, + arrayView1d< localIndex const > const bottomSurfaceNodeIndicator, + arrayView1d< real32 > referenceFrequencies, + arrayView1d< real32 > anelasticityCoefficients ) + + { + real64 const dt2 = pow( dt, 2 ); + integer nl = referenceFrequencies.size( 0 ); + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + if( freeSurfaceNodeIndicator[a] != 1 ) + { + p_np1[a] = 2.0*mass[a]*p_n[a]/dt2; + p_np1[a] -= mass[a]*p_nm1[a]/dt2; + p_np1[a] += stiffnessVector_p[a]; + p_np1[a] += rhs[a]; + + q_np1[a] = 2.0*mass[a]*q_n[a]/dt2; + q_np1[a] -= mass[a]*q_nm1[a]/dt2; + q_np1[a] += stiffnessVector_q[a]; + q_np1[a] += rhs[a]; + + // Apply attenuation + for( integer l = 0; l < nl; l++ ) + { + real32 gammal = ( 2.0 - referenceFrequencies[ l ] * dt )/( 2.0 + referenceFrequencies[ l ] * dt ); + real32 betal = anelasticityCoefficients[ l ] * referenceFrequencies[ l ] * 2.0 * dt /( 2.0 + referenceFrequencies[ l ] * dt ); + real32 gammalp = 0.5 + 0.5 * gammal; + real32 betalp = 0.5 * betal; + p_np1[a] += ( gammalp * divpsi_p( a, l ) + betalp * stiffnessVectorA_p[ a ] ); + q_np1[a] += ( gammalp * divpsi_q( a, l ) + betalp * stiffnessVectorA_q[ a ] ); + divpsi_p( a, l ) = gammal * divpsi_p( a, l ) + betal * stiffnessVectorA_p[ a ]; + divpsi_q( a, l ) = gammal * divpsi_q( a, l ) + betal * stiffnessVectorA_q[ a ]; + } + + if( lateralSurfaceNodeIndicator[a] != 1 && bottomSurfaceNodeIndicator[a] != 1 ) + { + // Interior node, no boundary terms + p_np1[a] /= mass[a]/dt2; + q_np1[a] /= mass[a]/dt2; + } + else + { + // Boundary node + p_np1[a] += damping_p[a]*p_nm1[a]/dt/2; + p_np1[a] += damping_pq[a]*q_nm1[a]/dt/2; + + q_np1[a] += damping_q[a]*q_nm1[a]/dt/2; + q_np1[a] += damping_qp[a]*p_nm1[a]/dt/2; + // Hand-made Inversion of 2x2 matrix + real32 coef_pp = mass[a]/dt2; + coef_pp += damping_p[a]/dt/2; + real32 coef_pq = damping_pq[a]/dt/2; + + real32 coef_qq = mass[a]/dt2; + coef_qq += damping_q[a]/2/dt; + real32 coef_qp = damping_qp[a]/dt/2; + + real32 det_pq = 1/(coef_pp * coef_qq - coef_pq*coef_qp); + + real32 aux_p_np1 = p_np1[a]; + p_np1[a] = det_pq*(coef_qq*p_np1[a] - coef_pq*q_np1[a]); + q_np1[a] = det_pq*(coef_pp*q_np1[a] - coef_qp*aux_p_np1); + } + } + } ); + }; }; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp index acf9572b21..2fa8ac11a6 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp @@ -1,4 +1,4 @@ -/* + /* * ------------------------------------------------------------------------------------------------------------ * SPDX-License-Identifier: LGPL-2.1-only * @@ -439,26 +439,13 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() dampingz ); } ); } ); - - // check anelasticity coefficient and/or compute it if needed - if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) - { - real32 minQVal = computeGlobalMinQFactor(); - if( m_slsAnelasticityCoefficients.size( 0 ) == 1 && m_slsAnelasticityCoefficients[ 0 ] < 0 ) - { - m_slsAnelasticityCoefficients[ 0 ] = 2.0 * minQVal / ( minQVal - 1.0 ); - } - // test if anelasticity is too high and artifacts could appear - real32 ySum = 0.0; - for( integer l = 0; l < m_slsAnelasticityCoefficients.size( 0 ); l++ ) - { - ySum += m_slsAnelasticityCoefficients[ l ]; - } - GEOS_WARNING_IF( ySum > minQVal, "The anelasticity parameters are too high for the given quality factor. This could lead to solution artifacts such as zero-velocity waves." ); - } - } ); + // check anelasticity coefficient and/or compute it if needed + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + initializeAnelasticityCoefficients< elasticfields::ElasticQualityFactorP, elasticfields::ElasticQualityFactorS >( domain ); + } if( m_timestepStabilityLimit==1 ) { real64 dtOut = 0.0; @@ -631,30 +618,6 @@ real64 ElasticWaveEquationSEM::computeTimeStep( real64 & dtOut ) return m_timeStep * m_cflFactor; } -real32 ElasticWaveEquationSEM::computeGlobalMinQFactor() -{ - RAJA::ReduceMin< ReducePolicy< EXEC_POLICY >, real32 > minQ( LvArray::NumericLimits< real32 >::max ); - DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & elementSubRegion ) - { - arrayView1d< real32 const > const qp = elementSubRegion.getField< elasticfields::ElasticQualityFactorP >(); - arrayView1d< real32 const > const qs = elementSubRegion.getField< elasticfields::ElasticQualityFactorS >(); - forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const e ) { - minQ.min( qp[e] ); - minQ.min( qs[e] ); - } ); - } ); - } ); - real32 minQVal = minQ.get(); - return MpiWrapper::min< real32 >( minQVal ); -} - void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartition & domain ) { FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp index 8907f35cd7..b0487136ed 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp @@ -157,12 +157,6 @@ class ElasticWaveEquationSEM : public WaveSolverBase void prepareNextTimestep( MeshLevel & mesh ); - /** - * @brief Computes the minimum attenuation quality factor over all the mesh. This is useful for computing anelasticity coefficients, which - * are usually global parameters - */ - real32 computeGlobalMinQFactor(); - protected: virtual void postInputInitialization() override final; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/shared/ElasticTimeSchemeSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/shared/ElasticTimeSchemeSEMKernel.hpp index 9b36fea044..f39281867a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/shared/ElasticTimeSchemeSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/shared/ElasticTimeSchemeSEMKernel.hpp @@ -128,6 +128,8 @@ struct ElasticTimeSchemeSEM * @param[in] rhsy the right-hand-side for displacement in y-direction * @param[in] rhsz the right-hand-side for displacement in z-direction * @param[in] solverTargetNodesSet the targetted nodeset (useful in particular when we do elasto-acoustic simulation ) + * @param[in] referenceFrequencies the reference frequencies for each SLS + * @param[in] anelasticityCoefficients the coefficients of anelasticity for each SLS */ static void AttenuationLeapFrog( real64 const dt, arrayView1d< real32 > const ux_np1, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp index b55ca1f38d..122048df11 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp @@ -27,8 +27,6 @@ #include "fieldSpecification/PerfectlyMatchedLayer.hpp" #include "mainInterface/ProblemManager.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "mesh/DomainPartition.hpp" -#include "WaveSolverUtils.hpp" #include "events/EventManager.hpp" #include @@ -565,5 +563,4 @@ bool WaveSolverBase::directoryExists( std::string const & directoryName ) return stat( directoryName.c_str(), &buffer ) == 0; } - } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp index 2c319ca4f3..091b70f810 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp @@ -29,6 +29,7 @@ #if !defined( GEOS_USE_HIP ) #include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" #endif +#include "mesh/DomainPartition.hpp" #include "WaveSolverUtils.hpp" #if !defined( GEOS_USE_HIP ) @@ -143,6 +144,58 @@ class WaveSolverBase : public PhysicsSolverBase void computeTargetNodeSet( arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes, localIndex const subRegionSize, localIndex const numQuadraturePointsPerElem ); + + /** + * @brief Computes the minimum value of the given element field over the whole mesh. + * @param[in] F The type of the field whose minimum must be computed + * @param[in] T The data type of the field + * @eturn the minimum value + */ + template< typename F, typename T > + T computeGlobalMinOfElemField( DomainPartition & domain ) + { + RAJA::ReduceMin< ReducePolicy< EXEC_POLICY >, T > minF( LvArray::NumericLimits< T >::max ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< T const > const f = elementSubRegion.getField< F >(); + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const e ) { + minF.min( f[e] ); + } ); + } ); + } ); + T minFVal = minF.get(); + return MpiWrapper::min< T >( minFVal ); + } + + /** + * @brief Initializes the anealsticity coefficients if needed, and checks that the anealsticity values + * are within a valid range. Gives a warning if not. + * @param[in] Qs the quality factor fields used to compute the anelasticity value if not given + */ + template< typename... Qs > + void initializeAnelasticityCoefficients( DomainPartition & domain ) + { + // compute miniumum quality factor + real32 minQVal = LvArray::NumericLimits< real32>::max; + ( (minQVal = std::min( minQVal, computeGlobalMinOfElemField< Qs, real32 >( domain ) ) ), ... ); + if( m_slsAnelasticityCoefficients.size( 0 ) == 1 && m_slsAnelasticityCoefficients[ 0 ] < 0 ) + { + m_slsAnelasticityCoefficients[ 0 ] = 2.0 * minQVal / ( minQVal - 1.0 ); + } + // test if anelasticity is too high and artifacts could appear + real32 ySum = 0.0; + for( integer l = 0; l < m_slsAnelasticityCoefficients.size( 0 ); l++ ) + { + ySum += m_slsAnelasticityCoefficients[ l ]; + } + GEOS_WARNING_IF( ySum > minQVal, "The anelasticity parameters are too high for the given quality factor. This could lead to solution artifacts such as zero-velocity waves." ); + } protected: diff --git a/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt b/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt index c1296f06b1..9b0bbac73c 100644 --- a/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt @@ -5,7 +5,9 @@ set( gtest_geosx_tests testWavePropagationElasticFirstOrder.cpp testWavePropagationDAS.cpp testWavePropagationElasticVTI.cpp - testWavePropagationAttenuation.cpp + testWavePropagationAttenuationElastic.cpp + testWavePropagationAttenuationAcoustic.cpp + testWavePropagationAttenuationAcousticVTI.cpp testWavePropagationAcousticFirstOrder.cpp testWavePropagationAdjoint1.cpp ) diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcoustic.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcoustic.cpp new file mode 100644 index 0000000000..fee7280d38 --- /dev/null +++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcoustic.cpp @@ -0,0 +1,222 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// using some utility classes from the following unit test +#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" + +#include "common/DataTypes.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/DomainPartition.hpp" +#include "mainInterface/GeosxState.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp" +#include "physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp" + +#include + +using namespace geos; +using namespace geos::dataRepository; +using namespace geos::testing; + +CommandLineOptions g_commandLineOptions; + +// This unit test checks the interpolation done to extract seismic traces from a wavefield. +// It computes a seismogram at a receiver co-located with the source and compares it to the surrounding receivers. +char const * xmlInput = + R"xml( + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + )xml"; + +class AcousticWaveEquationSEMTest : public ::testing::Test +{ +public: + + AcousticWaveEquationSEMTest(): + state( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ) + {} + +protected: + + void SetUp() override + { + setupProblemFromXML( state.getProblemManager(), xmlInput ); + } + + static real64 constexpr time = 0.0; + static real64 constexpr dt = 1e-1; + static real64 constexpr eps = std::numeric_limits< real64 >::epsilon(); + + GeosxState state; + AcousticWaveEquationSEM * propagator; +}; + +real64 constexpr AcousticWaveEquationSEMTest::time; +real64 constexpr AcousticWaveEquationSEMTest::dt; +real64 constexpr AcousticWaveEquationSEMTest::eps; + +TEST_F( AcousticWaveEquationSEMTest, SeismoTrace ) +{ + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + propagator = &state.getProblemManager().getPhysicsSolverManager().getGroup< AcousticWaveEquationSEM >( "acousticSolver" ); + real64 time_n = time; + // run for 1s (10 steps) + for( int i=0; i<10; i++ ) + { + propagator->explicitStepForward( time_n, dt, i, domain, false ); + time_n += dt; + } + // cleanup (triggers calculation of the remaining seismograms data points) + propagator->cleanup( 1.0, 10, 0, 0, domain ); + + // retrieve seismo + arrayView2d< real32 > const pressureReceivers = propagator->getReference< array2d< real32 > >( AcousticWaveEquationSEM::viewKeyStruct::pressureNp1AtReceiversString() ).toView(); + + // move it to CPU, if needed + pressureReceivers.move( hostMemorySpace, false ); + + // check number of seismos and trace length + ASSERT_EQ( pressureReceivers.size( 1 ), 10 ); + ASSERT_EQ( pressureReceivers.size( 0 ), 11 ); + + // check pressure content. The signal values cannot be directly checked as the problem is too small. + // Since the basis is linear, check that the seismograms are nonzero (for t>0) and the seismogram at the center is equal + // to the average of the others. + for( int i = 0; i < 11; i++ ) + { + if( i > 0 ) + { + ASSERT_TRUE( std::abs( pressureReceivers[i][8] ) > 0 ); + } + double avg = 0; + for( int r=0; r<8; r++ ) + { + avg += pressureReceivers[i][r]; + } + avg /= 8.0; + ASSERT_TRUE( std::abs( pressureReceivers[i][8] - avg ) < 0.00001 ); + } +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; +} diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcousticVTI.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcousticVTI.cpp new file mode 100644 index 0000000000..090f5f4f71 --- /dev/null +++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcousticVTI.cpp @@ -0,0 +1,236 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// using some utility classes from the following unit test +#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" + +#include "common/DataTypes.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/DomainPartition.hpp" +#include "mainInterface/GeosxState.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp" +#include "physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp" + +#include + +using namespace geos; +using namespace geos::dataRepository; +using namespace geos::testing; + +CommandLineOptions g_commandLineOptions; + +// This unit test checks the interpolation done to extract seismic traces from a wavefield. +// It computes a seismogram at a receiver co-located with the source and compares it to the surrounding receivers. +char const * xmlInput = + R"xml( + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + )xml"; + +class AcousticVTIWaveEquationSEMTest : public ::testing::Test +{ +public: + + AcousticVTIWaveEquationSEMTest(): + state( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ) + {} + +protected: + + void SetUp() override + { + setupProblemFromXML( state.getProblemManager(), xmlInput ); + } + + static real64 constexpr time = 0.0; + static real64 constexpr dt = 1e-1; + static real64 constexpr eps = std::numeric_limits< real64 >::epsilon(); + + GeosxState state; + AcousticVTIWaveEquationSEM * propagator; +}; + +real64 constexpr AcousticVTIWaveEquationSEMTest::time; +real64 constexpr AcousticVTIWaveEquationSEMTest::dt; +real64 constexpr AcousticVTIWaveEquationSEMTest::eps; + +TEST_F( AcousticVTIWaveEquationSEMTest, SeismoTrace ) +{ + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + propagator = &state.getProblemManager().getPhysicsSolverManager().getGroup< AcousticVTIWaveEquationSEM >( "acousticVTISolver" ); + real64 time_n = time; + // run for 1s (10 steps) + for( int i=0; i<10; i++ ) + { + propagator->explicitStepForward( time_n, dt, i, domain, false ); + time_n += dt; + } + // cleanup (triggers calculation of the remaining seismograms data points) + propagator->cleanup( 1.0, 10, 0, 0, domain ); + + // retrieve seismo + arrayView2d< real32 > const pressureReceivers = propagator->getReference< array2d< real32 > >( AcousticVTIWaveEquationSEM::viewKeyStruct::pressureNp1AtReceiversString() ).toView(); + + // move it to CPU, if needed + pressureReceivers.move( hostMemorySpace, false ); + + // check number of seismos and trace length + ASSERT_EQ( pressureReceivers.size( 1 ), 10 ); + ASSERT_EQ( pressureReceivers.size( 0 ), 11 ); + + // check pressure content. The signal values cannot be directly checked as the problem is too small. + // Since the basis is linear, check that the seismograms are nonzero (for t>0) and the seismogram at the center is equal + // to the average of the others. + for( int i = 0; i < 11; i++ ) + { + if( i > 0 ) + { + ASSERT_TRUE( std::abs( pressureReceivers[i][8] ) > 0 ); + } + double avg = 0; + for( int r=0; r<8; r++ ) + { + avg += pressureReceivers[i][r]; + } + avg /= 8.0; + ASSERT_TRUE( std::abs( pressureReceivers[i][8] - avg ) < 0.00001 ); + } +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; +} diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuation.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationElastic.cpp similarity index 100% rename from src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuation.cpp rename to src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationElastic.cpp