From 065a9a1eb9444a149415f6e9e3335df59ab3df61 Mon Sep 17 00:00:00 2001 From: j0405284 Date: Mon, 14 Oct 2024 12:26:12 +0200 Subject: [PATCH 1/7] first version of viscoacoustic, isotropic and VTI (no PML) --- .../AcousticVTIWaveEquationSEM.cpp | 43 ++++- .../isotropic/AcousticWaveEquationSEM.cpp | 60 ++++++- .../AcousticWaveEquationSEMKernel.hpp | 107 ++++++++++-- .../sem/acoustic/shared/AcousticFields.hpp | 56 +++++++ .../shared/AcousticTimeSchemeSEMKernel.hpp | 156 ++++++++++++++++++ 5 files changed, 397 insertions(+), 25 deletions(-) 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 edf012c68d8..c10423fa86c 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< elasticfields::DivPsi_p, + elasticfields::DivPsi_q, + elasticfields::StiffnessVectorA_p, + elasticfields::StiffnessVectorA_q >( getName() ); + nodeManager.getField< elasticfields::DivPsi_p >().resizeDimension< 1 >( l ); + nodeManager.getField< elasticfields::DivPsi_q >().resizeDimension< 1 >( l ); + } FaceManager & faceManager = mesh.getFaceManager(); faceManager.registerField< acousticfields::AcousticFreeSurfaceFaceIndicator >( getName() ); @@ -608,11 +618,25 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, 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< elasticfields::StiffnessVectorA_p >(); + arrayView1d< real32 > const stiffnessVectorA_q = nodeManager.getField< elasticfields::StiffnessVectorA_q >(); + arrayView2d< real32 > const divpsix = nodeManager.getField< elasticfields::DivPsi_p >(); + arrayView2d< real32 > const divpsiy = nodeManager.getField< elasticfields::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, mass, stiffnessVector_p, + stiffnessVector_q, stiffnessVectorA_p, stiffnessVectorA_q, damping_p, damping_pq, damping_q, + damping_qp, rhs, freeSurfaceNodeIndicator, lateralSurfaceNodeIndicator, bottomSurfaceNodeIndicator ); + } + 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() } ); @@ -637,6 +661,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< elasticfields::StiffnessVectorA_p >(); + arrayView1d< real32 > const stiffnessVectorA_q = nodeManager.getField< elasticfields::StiffnessVectorA_q >(); + forAll< EXEC_POLICY >( nodeManager.size() , [=] GEOS_HOST_DEVICE ( localIndex const q ) + { + stiffnessVectorA_p[a] = stiffnessVectorA_q[a] = 0.0; + } ); + } } ); 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 6e1852d0a55..9aabd85b02a 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< elasticfields::DivPsi, + elasticfields::StiffnessVectorA >( getName() ); + nodeManager.getField< elasticfields::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() ); + } } ); } ); @@ -337,10 +349,25 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() velocity, density, damping ); - - } ); } ); + + // 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." ); + } } ); WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); @@ -943,6 +970,15 @@ void AcousticWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh ) stiffnessVector[a] = rhs[a] = 0.0; } ); + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< elasticfields::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, @@ -985,11 +1021,27 @@ 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 ) + { + AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, + rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); + } + else + { + arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< elasticfields::StiffnessVectorA >(); + arrayView2d< real32 > const divpsi = nodeManager.getField< elasticfields::DivPsi >(); + arrayView1d< real32 > const referenceFrequencies = m_slsReferenceAngularFrequencies.toView(); + arrayView1d< real32 > const anelasticityCoefficients = m_slsAnelasticityCoefficients.toView(); + AcousticTimeSchemeSEM::AttenuationLeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, stiffnessVectorA, + damping, rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); + } } else { 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 f582a027f97..e10f867c8f5 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 ) { @@ -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::StiffnessVector >; + +//***************************************************************************** + /** + * @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::AcusticQualityFactor >() ), + {} + + /** + * @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 dd067a11a1f..47557db5339 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( PressureDoubleDerivative, "pressureDoubleDerivative", array1d< real32 >, @@ -123,6 +147,30 @@ 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( 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( DampingVector, "dampingVector", array1d< real32 >, @@ -147,6 +195,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 20b838f6d62..542f4ecd70e 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp @@ -69,6 +69,57 @@ 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] 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 ) + */ + static void AttenuationLeapFrogWithoutPML( real64 const dt, + arrayView1d< real32 > const p_np1, + arrayView1d< real32 > const p_n, + arrayView1d< real32 > const p_nm1, + 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 ) + { + real64 const dt2 = pow( dt, 2 ); + 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 @@ -127,6 +178,19 @@ struct AcousticTimeSchemeSEM 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( a, l ) + betalp * stiffnessVectorA_p[ a ] ); + q_np1[a] += ( gammalp * divpsi( 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 @@ -160,6 +224,98 @@ 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] 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 + */ + static void LeapFrogforVTI( 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, + arrayView1d< real32 const > const mass, + arrayView1d< real32 > const stiffnessVector_p, + arrayView1d< real32 > const stiffnessVector_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 ) + + { + real64 const dt2 = pow( dt, 2 ); + 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]; + + 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); + } + } + } ); + }; }; From dfb12a56fe1a6e31847aa8e926bd8df0cd023150 Mon Sep 17 00:00:00 2001 From: j0405284 Date: Tue, 15 Oct 2024 17:07:02 +0200 Subject: [PATCH 2/7] addd a lot of leftovers --- .../anisotropic/AcousticVTIFields.hpp | 32 +++++ .../AcousticVTIWaveEquationSEM.cpp | 53 +++++--- .../AcousticVTIWaveEquationSEMKernel.hpp | 125 ++++++++++++++---- .../isotropic/AcousticWaveEquationSEM.cpp | 52 ++++---- .../AcousticWaveEquationSEMKernel.hpp | 10 +- .../sem/acoustic/shared/AcousticFields.hpp | 32 ----- .../shared/AcousticTimeSchemeSEMKernel.hpp | 84 +++++++----- .../isotropic/ElasticWaveEquationSEM.cpp | 47 +------ .../isotropic/ElasticWaveEquationSEM.hpp | 6 - .../shared/ElasticTimeSchemeSEMKernel.hpp | 2 + .../wavePropagation/shared/WaveSolverBase.cpp | 44 ++++++ .../wavePropagation/shared/WaveSolverBase.hpp | 18 +++ 12 files changed, 323 insertions(+), 182 deletions(-) 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 1069ee8ae31..5e26f6b25e6 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 c10423fa86c..1ad6668d3b6 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp @@ -87,12 +87,12 @@ void AcousticVTIWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { integer l = m_slsReferenceAngularFrequencies.size( 0 ); - nodeManager.registerField< elasticfields::DivPsi_p, - elasticfields::DivPsi_q, - elasticfields::StiffnessVectorA_p, - elasticfields::StiffnessVectorA_q >( getName() ); - nodeManager.getField< elasticfields::DivPsi_p >().resizeDimension< 1 >( l ); - nodeManager.getField< elasticfields::DivPsi_q >().resizeDimension< 1 >( l ); + 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(); @@ -341,6 +341,12 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); + // check anelasticity coefficient and/or compute it if needed + if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) + { + initializeAnelasticityCoefficients< acousticfields::AcousticQualityFactor >(); + } + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); } @@ -612,6 +618,20 @@ 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( cycleNumber, rhs ); /// calculate your time integrators @@ -620,15 +640,16 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { - arrayView1d< real32 > const stiffnessVectorA_p = nodeManager.getField< elasticfields::StiffnessVectorA_p >(); - arrayView1d< real32 > const stiffnessVectorA_q = nodeManager.getField< elasticfields::StiffnessVectorA_q >(); - arrayView2d< real32 > const divpsix = nodeManager.getField< elasticfields::DivPsi_p >(); - arrayView2d< real32 > const divpsiy = nodeManager.getField< elasticfields::DivPsi_q >(); + 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, mass, stiffnessVector_p, - stiffnessVector_q, stiffnessVectorA_p, stiffnessVectorA_q, damping_p, damping_pq, damping_q, - damping_qp, rhs, freeSurfaceNodeIndicator, lateralSurfaceNodeIndicator, bottomSurfaceNodeIndicator ); + 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 { @@ -663,9 +684,9 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, } ); if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { - arrayView1d< real32 > const stiffnessVectorA_p = nodeManager.getField< elasticfields::StiffnessVectorA_p >(); - arrayView1d< real32 > const stiffnessVectorA_q = nodeManager.getField< elasticfields::StiffnessVectorA_q >(); - forAll< EXEC_POLICY >( nodeManager.size() , [=] GEOS_HOST_DEVICE ( localIndex const q ) + 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 3a711540091..b96545b986e 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp @@ -390,12 +390,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: @@ -429,22 +430,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 >() ), @@ -459,7 +460,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 @@ -475,6 +476,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 ]{}; @@ -501,6 +503,7 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE, stack.xLocal[ a ][ i ] = m_nodeCoords[ nodeIndex ][ i ]; } } + stack.factor = 1.0; } GEOS_HOST_DEVICE @@ -518,7 +521,7 @@ class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE, /** * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel * - * ### ExplicitAcousticVTISEM Description + * ### ExplicitAcousticVTISEMBase Description * Calculates stiffness vector * */ @@ -532,9 +535,9 @@ 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 @@ -542,10 +545,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; } ); } @@ -581,10 +584,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 9aabd85b02a..c8b6f73b137 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -87,9 +87,9 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { integer l = m_slsReferenceAngularFrequencies.size( 0 ); - nodeManager.registerField< elasticfields::DivPsi, - elasticfields::StiffnessVectorA >( getName() ); - nodeManager.getField< elasticfields::DivPsi >().resizeDimension< 1 >( l ); + 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 @@ -351,25 +351,14 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() damping ); } ); } ); - - // 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< acousticfields::AcousticQualityFactor >(); + } + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); } @@ -972,7 +961,7 @@ void AcousticWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh ) } ); if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { - arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< elasticfields::StiffnessVectorA >(); + arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< acousticfields::StiffnessVectorA >(); forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n ) { localIndex const a = solverTargetNodesSet[n]; @@ -1011,6 +1000,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 EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() ); @@ -1035,12 +1036,13 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, } else { - arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< elasticfields::StiffnessVectorA >(); - arrayView2d< real32 > const divpsi = nodeManager.getField< elasticfields::DivPsi >(); + 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, mass, stiffnessVector, stiffnessVectorA, - damping, rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); + AcousticTimeSchemeSEM::AttenuationLeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, divpsi, mass, stiffnessVector, stiffnessVectorA, + damping, rhs, freeSurfaceNodeIndicator, solverTargetNodesSet, referenceFrequencies, + anelasticityCoefficients ); } } else 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 e10f867c8f5..2e7c9f761fe 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEMKernel.hpp @@ -119,7 +119,7 @@ class ExplicitAcousticSEMBase : 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 @@ -179,7 +179,7 @@ class ExplicitAcousticSEMBase : public finiteElement::KernelBase< SUBREGION_TYPE /** * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel * - * ### ExplicitAcousticSEM Description + * ### ExplicitAcousticSEMBase Description * Calculates stiffness vector * */ @@ -241,7 +241,7 @@ class ExplicitAcousticAttenuativeSEM : public ExplicitAcousticSEMBase< SUBREGION using Base = ExplicitAcousticSEMBase< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE, - fields::acousticfields::StiffnessVector >; + fields::acousticfields::StiffnessVectorA >; //***************************************************************************** /** @@ -269,7 +269,7 @@ class ExplicitAcousticAttenuativeSEM : public ExplicitAcousticSEMBase< SUBREGION finiteElementSpace, inputConstitutiveType, dt ), - m_qualityFactor( elementSubRegion.template getField< fields::acousticfields::AcusticQualityFactor >() ), + m_qualityFactor( elementSubRegion.template getField< fields::acousticfields::AcousticQualityFactor >() ) {} /** @@ -283,7 +283,7 @@ class ExplicitAcousticAttenuativeSEM : public ExplicitAcousticSEMBase< SUBREGION typename Base::StackVariables & stack ) const { Base::setup( k, stack ); - stack.invDensity = stack.invDensity * m_qualityFactor[ k ]; + stack.invDensity = stack.invDensity / m_qualityFactor[ k ]; } protected: diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp index 47557db5339..fe3a9d8d59a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp @@ -67,22 +67,6 @@ DECLARE_FIELD( DivPsi, 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( PressureDoubleDerivative, "pressureDoubleDerivative", array1d< real32 >, @@ -155,22 +139,6 @@ DECLARE_FIELD( StiffnessVectorA, WRITE_AND_READ, "acoustic attenuation stiffness vector." ); -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( DampingVector, "dampingVector", array1d< real32 >, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp index 542f4ecd70e..23aa59a8b1d 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp @@ -75,6 +75,7 @@ struct AcousticTimeSchemeSEM * @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 @@ -82,20 +83,26 @@ struct AcousticTimeSchemeSEM * @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 ) + 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]; @@ -178,19 +185,6 @@ struct AcousticTimeSchemeSEM 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( a, l ) + betalp * stiffnessVectorA_p[ a ] ); - q_np1[a] += ( gammalp * divpsi( 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 @@ -234,6 +228,8 @@ struct AcousticTimeSchemeSEM * @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 @@ -247,29 +243,38 @@ struct AcousticTimeSchemeSEM * @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 LeapFrogforVTI( 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, - arrayView1d< real32 const > const mass, - arrayView1d< real32 > const stiffnessVector_p, - arrayView1d< real32 > const stiffnessVector_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 ) + 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 ) @@ -284,6 +289,19 @@ struct AcousticTimeSchemeSEM 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 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 b6e28c5f3d5..91be2e20ab3 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp @@ -454,53 +454,18 @@ 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 >(); + } + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); WaveSolverUtils::initTrace( "dasTraceReceiver", getName(), m_outputSeismoTrace, m_linearDASGeometry.size( 0 ), m_receiverIsLocal ); } -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 ) { 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 f3c0b14396f..903da7be76e 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 55f7dcc163d..49ef118b157 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 d15555e75a9..0371c4bb976 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp @@ -540,5 +540,49 @@ bool WaveSolverBase::directoryExists( std::string const & directoryName ) return stat( directoryName.c_str(), &buffer ) == 0; } +template< typename F, typename T > +real32 WaveSolverBase::computeGlobalMinOfElemField() +{ + RAJA::ReduceMin< ReducePolicy< EXEC_POLICY >, T > minF( LvArray::NumericLimits< T >::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< 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 ); +} + +template< typename... Qs > +void WaveSolverBase::initializeAnelasticityCoefficients() +{ + // compute miniumum quality factor + real32 minQVal = LvArray::NumericLimits< real32>::max; + ( minQVal = std::min( minQVal, computeGlobalMinOfElemField< Qs >() ), ... ); + 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." ); + +} + + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp index 6c6ce292243..32083c400a1 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp @@ -138,6 +138,24 @@ class WaveSolverBase : public SolverBase 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 > + real32 computeGlobalMinOfElemField(); + + /** + * @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(); + protected: From 2eaf425be6bedcd1b7108df61ec28a79096f269e Mon Sep 17 00:00:00 2001 From: j0405284 Date: Tue, 15 Oct 2024 18:01:13 +0200 Subject: [PATCH 3/7] missed a couple things --- .../AcousticVTIWaveEquationSEM.cpp | 2 +- .../isotropic/AcousticWaveEquationSEM.cpp | 2 +- .../isotropic/ElasticWaveEquationSEM.cpp | 2 +- .../wavePropagation/shared/WaveSolverBase.cpp | 47 ------------------- .../wavePropagation/shared/WaveSolverBase.hpp | 41 ++++++++++++++-- 5 files changed, 41 insertions(+), 53 deletions(-) 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 1ad6668d3b6..dcce3efd57b 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp @@ -344,7 +344,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups() // check anelasticity coefficient and/or compute it if needed if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { - initializeAnelasticityCoefficients< acousticfields::AcousticQualityFactor >(); + initializeAnelasticityCoefficients< acousticfields::AcousticQualityFactor >( domain ); } WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); 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 c8b6f73b137..2ba5baf5e1a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -356,7 +356,7 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() // check anelasticity coefficient and/or compute it if needed if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { - initializeAnelasticityCoefficients< acousticfields::AcousticQualityFactor >(); + initializeAnelasticityCoefficients< acousticfields::AcousticQualityFactor >( domain ); } WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); 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 91be2e20ab3..7c730087bee 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp @@ -459,7 +459,7 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() // check anelasticity coefficient and/or compute it if needed if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) { - initializeAnelasticityCoefficients< elasticfields::ElasticQualityFactorP, elasticfields::ElasticQualityFactorS >(); + initializeAnelasticityCoefficients< elasticfields::ElasticQualityFactorP, elasticfields::ElasticQualityFactorS >( domain ); } WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp index 0371c4bb976..b2aafc97039 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 @@ -540,49 +538,4 @@ bool WaveSolverBase::directoryExists( std::string const & directoryName ) return stat( directoryName.c_str(), &buffer ) == 0; } -template< typename F, typename T > -real32 WaveSolverBase::computeGlobalMinOfElemField() -{ - RAJA::ReduceMin< ReducePolicy< EXEC_POLICY >, T > minF( LvArray::NumericLimits< T >::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< 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 ); -} - -template< typename... Qs > -void WaveSolverBase::initializeAnelasticityCoefficients() -{ - // compute miniumum quality factor - real32 minQVal = LvArray::NumericLimits< real32>::max; - ( minQVal = std::min( minQVal, computeGlobalMinOfElemField< Qs >() ), ... ); - 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." ); - -} - - - } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp index 32083c400a1..156e2955e36 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp @@ -28,6 +28,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 ) @@ -146,7 +147,26 @@ class WaveSolverBase : public SolverBase * @eturn the minimum value */ template< typename F, typename T > - real32 computeGlobalMinOfElemField(); + 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 @@ -154,8 +174,23 @@ class WaveSolverBase : public SolverBase * @param[in] Qs the quality factor fields used to compute the anelasticity value if not given */ template< typename... Qs > - void initializeAnelasticityCoefficients(); - + 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: From 3becd25d458f42d1a96782cd6ff0ceab1b7595c5 Mon Sep 17 00:00:00 2001 From: j0405284 Date: Tue, 15 Oct 2024 19:25:36 +0200 Subject: [PATCH 4/7] fixed a few bugs --- .../anisotropic/AcousticVTIFields.hpp | 4 ++-- .../anisotropic/AcousticVTIWaveEquationSEM.cpp | 5 +++++ .../isotropic/AcousticWaveEquationSEM.cpp | 15 ++++++++++----- .../sem/acoustic/shared/AcousticFields.hpp | 2 +- 4 files changed, 18 insertions(+), 8 deletions(-) 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 5e26f6b25e6..791f3ed9029 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFields.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFields.hpp @@ -80,7 +80,7 @@ DECLARE_FIELD( StiffnessVectorA_p, 0, NOPLOT, WRITE_AND_READ, - "p-type acoustic attenuation stiffness vector." ); + "P-type acoustic attenuation stiffness vector." ); DECLARE_FIELD( StiffnessVectorA_q, "stiffnessVectorA_q", @@ -88,7 +88,7 @@ DECLARE_FIELD( StiffnessVectorA_q, 0, NOPLOT, WRITE_AND_READ, - "q-type acoustic attenuation stiffness vector." ); + "Q-type acoustic attenuation stiffness vector." ); DECLARE_FIELD( Pressure_p_nm1, "pressure_p_nm1", 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 dcce3efd57b..662450f06d4 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp @@ -663,6 +663,11 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, 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, 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 2ba5baf5e1a..f9522b9de45 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -1030,11 +1030,6 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, { GEOS_MARK_SCOPE ( updateP ); if( m_attenuationType == WaveSolverUtils::AttenuationType::sls ) - { - AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, - rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); - } - else { arrayView1d< real32 > const stiffnessVectorA = nodeManager.getField< acousticfields::StiffnessVectorA >(); arrayView2d< real32 > const divpsi = nodeManager.getField< acousticfields::DivPsi >(); @@ -1044,6 +1039,11 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, damping, rhs, freeSurfaceNodeIndicator, solverTargetNodesSet, referenceFrequencies, anelasticityCoefficients ); } + else + { + AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, + rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); + } } else { @@ -1128,6 +1128,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/shared/AcousticFields.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp index fe3a9d8d59a..e25364d7e7e 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp @@ -137,7 +137,7 @@ DECLARE_FIELD( StiffnessVectorA, 0, NOPLOT, WRITE_AND_READ, - "acoustic attenuation stiffness vector." ); + "Acoustic attenuation stiffness vector." ); DECLARE_FIELD( DampingVector, "dampingVector", From bc7d11fbfd1c76ca75fc8bbbf759f0b275ec503b Mon Sep 17 00:00:00 2001 From: j0405284 Date: Fri, 18 Oct 2024 14:34:51 +0200 Subject: [PATCH 5/7] added test and corrected bug --- .../AcousticVTIWaveEquationSEM.cpp | 4 + .../wavePropagationTests/CMakeLists.txt | 4 +- ...testWavePropagationAttenuationAcoustic.cpp | 222 ++++++++++++++++ ...tWavePropagationAttenuationAcousticVTI.cpp | 236 ++++++++++++++++++ ...testWavePropagationAttenuationElastic.cpp} | 0 5 files changed, 465 insertions(+), 1 deletion(-) create mode 100644 src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcoustic.cpp create mode 100644 src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcousticVTI.cpp rename src/coreComponents/unitTests/wavePropagationTests/{testWavePropagationAttenuation.cpp => testWavePropagationAttenuationElastic.cpp} (100%) 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 662450f06d4..90216ff3798 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp @@ -108,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() ); + } } ); } ); } diff --git a/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt b/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt index fab0a741e73..2adfe38f081 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 ) set( tplDependencyList ${parallelDeps} gtest ) diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcoustic.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAttenuationAcoustic.cpp new file mode 100644 index 00000000000..fee7280d38e --- /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 00000000000..090f5f4f719 --- /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 From dd43aa364c507dbaf9faefb5dab30a3daf287b63 Mon Sep 17 00:00:00 2001 From: j0405284 Date: Mon, 4 Nov 2024 09:37:45 +0100 Subject: [PATCH 6/7] leftover from merge --- .../secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp | 1 - 1 file changed, 1 deletion(-) 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 be3934d6ac8..f14a9daaf0b 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -352,7 +352,6 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ->>>>>>> develop WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); } From c9f1642d301fbb4cd607f2ec566a16a0166c537d Mon Sep 17 00:00:00 2001 From: j0405284 Date: Tue, 19 Nov 2024 16:15:35 +0100 Subject: [PATCH 7/7] leftover from merge --- .../isotropic/ElasticWaveEquationSEM.cpp | 27 +------------------ 1 file changed, 1 insertion(+), 26 deletions(-) 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 0b00c851bb6..368cf6a3870 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 * @@ -618,31 +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 ); -} ->>>>>>> develop - void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartition & domain ) { FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance();