diff --git a/src/coreComponents/physicsSolvers/wavePropagation/CMakeLists.txt b/src/coreComponents/physicsSolvers/wavePropagation/CMakeLists.txt index 54fe46da9fa..cdd45081727 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/wavePropagation/CMakeLists.txt @@ -4,6 +4,7 @@ set( physicsSolvers_headers wavePropagation/shared/WaveSolverBase.hpp wavePropagation/shared/WaveSolverUtils.hpp wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp + wavePropagation/shared/TaperKernel.hpp wavePropagation/sem/acoustic/shared/AcousticFields.hpp wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEMKernel.hpp @@ -27,7 +28,7 @@ set( physicsSolvers_headers wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcoustoElasticFields.hpp wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcoustoElasticTimeSchemeSEMKernel.hpp wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEM.hpp - wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEMKernel.hpp + wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEMKernel.hpp PARENT_SCOPE ) # Specify solver sources @@ -39,5 +40,5 @@ set( physicsSolvers_sources wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp - wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEM.cpp + wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEM.cpp PARENT_SCOPE) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp index 506c95a3d5d..431dcc4e07b 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp @@ -154,6 +154,28 @@ void AcousticFirstOrderWaveEquationSEM::postInputInitialization() m_receiverRegion.resize( numReceiversGlobal ); } +real32 AcousticFirstOrderWaveEquationSEM::getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + + real32 localMinWavespeed = 1e8; + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfields::AcousticVelocity >(); + real32 subRegionMinWavespeed = *std::min_element( velocity.begin(), velocity.end()); + if( localMinWavespeed > subRegionMinWavespeed ) + { + localMinWavespeed = subRegionMinWavespeed; + } + } ); + + real32 const globalMinWavespeed = MpiWrapper::min( localMinWavespeed ); + + return globalMinWavespeed; + +} + void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { GEOS_MARK_FUNCTION; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp index 173600e1b90..7fb78429816 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp @@ -69,6 +69,11 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase DomainPartition & domain, bool const computeGradient ) override; + /** + * @brief Get the minimum wavespeed on a mesh + */ + virtual real32 getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /** * @brief Initialize Perfectly Matched Layer (PML) information */ 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 28b2497012d..c7774b5f5b8 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp @@ -102,7 +102,6 @@ void AcousticVTIWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) } ); } - void AcousticVTIWaveEquationSEM::postInputInitialization() { @@ -113,6 +112,28 @@ void AcousticVTIWaveEquationSEM::postInputInitialization() m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); } +real32 AcousticVTIWaveEquationSEM::getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + + real32 localMinWavespeed = 1e8; + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfields::AcousticVelocity >(); + real32 subRegionMinWavespeed = *std::min_element( velocity.begin(), velocity.end()); + if( localMinWavespeed > subRegionMinWavespeed ) + { + localMinWavespeed = subRegionMinWavespeed; + } + } ); + + real32 const globalMinWavespeed = MpiWrapper::min( localMinWavespeed ); + + return globalMinWavespeed; + +} + void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp index dcfe6da4d14..88f98a4a161 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp @@ -69,6 +69,11 @@ class AcousticVTIWaveEquationSEM : public WaveSolverBase DomainPartition & GEOS_UNUSED_PARAM( domain ), bool const GEOS_UNUSED_PARAM( computeGradient ) ) override; + /** + * @brief Get the minimum wavespeed on a mesh + */ + virtual real32 getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /**@}*/ /** 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 c47e6e83fc1..23c810f75dc 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -33,6 +33,7 @@ #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp" #include "events/EventManager.hpp" #include "AcousticPMLSEMKernel.hpp" +#include "physicsSolvers/wavePropagation/shared/TaperKernel.hpp" #include "physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp" namespace geos @@ -120,9 +121,32 @@ void AcousticWaveEquationSEM::postInputInitialization() m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, m_receiverCoordinates.size( 0 ) + 1 ); } +real32 AcousticWaveEquationSEM::getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + + real32 localMinWavespeed = 1e8; + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfields::AcousticVelocity >(); + real32 subRegionMinWavespeed = *std::min_element( velocity.begin(), velocity.end()); + if( localMinWavespeed > subRegionMinWavespeed ) + { + localMinWavespeed = subRegionMinWavespeed; + } + } ); + + real32 const globalMinWavespeed = MpiWrapper::min( localMinWavespeed ); + + return globalMinWavespeed; + +} + void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { + GEOS_MARK_FUNCTION; arrayView1d< globalIndex const > const nodeLocalToGlobalMap = baseMesh.getNodeManager().localToGlobalMap().toViewConst(); @@ -328,15 +352,38 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); + if( m_timestepStabilityLimit==1 ) + { + real64 dtOut = 0.0; + computeTimeStep( dtOut ); + m_timestepStabilityLimit = 0; + m_timeStep=dtOut; + } + else + { + EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); + for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent ) + { + EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] ); + if( subEvent->getEventName() == "/Solvers/" + getName() ) + { + m_timeStep = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() ); + } + } + + } + + if( m_useTaper==1 ) + { + real32 vMin; + vMin = getGlobalMinWavespeed( mesh, regionNames ); + + arrayView1d< real32 > const taperCoeff = nodeManager.getField< fields::taperCoeff >(); + TaperKernel::computeTaperCoeff< EXEC_POLICY >( nodeManager.size(), nodeCoords, m_thicknessTaper, m_timeStep, vMin, m_reflectivityCoeff, + taperCoeff ); + } } ); - if( m_timestepStabilityLimit==1 ) - { - real64 dtOut = 0.0; - computeTimeStep( dtOut ); - m_timestepStabilityLimit = 0; - m_timeStep=dtOut; - } WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); @@ -1089,6 +1136,8 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, arrayView1d< real32 > const p_n = nodeManager.getField< acousticfields::Pressure_n >(); arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticfields::Pressure_np1 >(); + arrayView1d< real32 > const taperCoeff = nodeManager.getField< fields::taperCoeff >(); + arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfields::AcousticFreeSurfaceNodeIndicator >(); arrayView1d< real32 > const stiffnessVector = nodeManager.getField< acousticfields::StiffnessVector >(); arrayView1d< real32 > const rhs = nodeManager.getField< acousticfields::ForcingRHS >(); @@ -1115,6 +1164,13 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n, GEOS_MARK_SCOPE ( updateP ); AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); + + + if( m_useTaper==1 ) + { + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, p_np1 ); + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, p_n ); + } } else { diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp index 5f3f1c70885..365832402f5 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp @@ -80,6 +80,11 @@ class AcousticWaveEquationSEM : public WaveSolverBase DomainPartition & domain, bool const computeGradient ) override; + /** + * @brief Get the minimum wavespeed on a mesh + */ + virtual real32 getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /**@}*/ /** diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp index 5222bc115a3..d957f859ce2 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp @@ -202,7 +202,6 @@ DECLARE_FIELD( AuxiliaryVar4PML, NOPLOT, WRITE_AND_READ, "PML scalar auxiliary variable 4." ); - } } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp index 5569e40dcdf..608e2e337d4 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp @@ -202,6 +202,27 @@ void ElasticFirstOrderWaveEquationSEM::postInputInitialization() m_sigmayzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); } +real32 ElasticFirstOrderWaveEquationSEM::getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + + real32 localMinWavespeed = 1e8; + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real32 const > const velocity = elementSubRegion.getField< elasticfields::ElasticVelocityVs >(); + real32 subRegionMinWavespeed = *std::min_element( velocity.begin(), velocity.end()); + if( localMinWavespeed > subRegionMinWavespeed ) + { + localMinWavespeed = subRegionMinWavespeed; + } + } ); + + real32 const globalMinWavespeed = MpiWrapper::min( localMinWavespeed ); + + return globalMinWavespeed; + +} void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp index 4cdf886919a..a3daab96265 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp @@ -70,6 +70,11 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase DomainPartition & domain, bool const computeGradient ) override; + /** + * @brief Get the minimum wavespeed on a mesh (S-wavespeed in the elastic case) + */ + virtual real32 getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /** * @brief Initialize Perfectly Matched Layer (PML) information 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 acf9572b216..3a26e67681d 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp @@ -32,6 +32,7 @@ #include "physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp" #include "physicsSolvers/wavePropagation/sem/elastic/shared/ElasticTimeSchemeSEMKernel.hpp" #include "physicsSolvers/wavePropagation/sem/elastic/shared/ElasticMatricesSEMKernel.hpp" +#include "physicsSolvers/wavePropagation/shared/TaperKernel.hpp" #include "events/EventManager.hpp" namespace geos @@ -223,6 +224,27 @@ void ElasticWaveEquationSEM::postInputInitialization() } } +real32 ElasticWaveEquationSEM::getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + + real32 localMinWavespeed = 1e8; + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real32 const > const velocity = elementSubRegion.getField< elasticfields::ElasticVelocityVs >(); + real32 subRegionMinWavespeed = *std::min_element( velocity.begin(), velocity.end()); + if( localMinWavespeed > subRegionMinWavespeed ) + { + localMinWavespeed = subRegionMinWavespeed; + } + } ); + + real32 const globalMinWavespeed = MpiWrapper::min( localMinWavespeed ); + + return globalMinWavespeed; + +} void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { @@ -457,15 +479,40 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() 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." ); } + if( m_timestepStabilityLimit==1 ) + { + real64 dtOut = 0.0; + computeTimeStep( dtOut ); + m_timestepStabilityLimit = 0; + m_timeStep=dtOut; + } + else + { + EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); + for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent ) + { + EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] ); + if( subEvent->getEventName() == "/Solvers/" + getName() ) + { + m_timeStep = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() ); + } + } + + } + + if( m_useTaper==1 ) + { + real32 vMin; + vMin = getGlobalMinWavespeed( mesh, regionNames ); + + arrayView1d< real32 > const taperCoeff = nodeManager.getField< fields::taperCoeff >(); + TaperKernel::computeTaperCoeff< EXEC_POLICY >( nodeManager.size(), nodeCoords, m_thicknessTaper, m_timeStep, vMin, m_reflectivityCoeff, + taperCoeff ); + } + } ); - if( m_timestepStabilityLimit==1 ) - { - real64 dtOut = 0.0; - computeTimeStep( dtOut ); - m_timestepStabilityLimit = 0; - m_timeStep=dtOut; - } + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); WaveSolverUtils::initTrace( "dasTraceReceiver", getName(), m_outputSeismoTrace, m_linearDASGeometry.size( 0 ), m_receiverIsLocal ); @@ -778,6 +825,8 @@ void ElasticWaveEquationSEM::computeUnknowns( real64 const & time_n, arrayView1d< real32 > const uy_np1 = nodeManager.getField< elasticfields::Displacementy_np1 >(); arrayView1d< real32 > const uz_np1 = nodeManager.getField< elasticfields::Displacementz_np1 >(); + arrayView1d< real32 > const taperCoeff = nodeManager.getField< fields::taperCoeff >(); + arrayView1d< real32 > const rhsx = nodeManager.getField< elasticfields::ForcingRHSx >(); arrayView1d< real32 > const rhsy = nodeManager.getField< elasticfields::ForcingRHSy >(); arrayView1d< real32 > const rhsz = nodeManager.getField< elasticfields::ForcingRHSz >(); @@ -846,6 +895,18 @@ void ElasticWaveEquationSEM::computeUnknowns( real64 const & time_n, mass, dampingx, dampingy, dampingz, stiffnessVectorx, stiffnessVectory, stiffnessVectorz, rhsx, rhsy, rhsz, solverTargetNodesSet ); } + + if( m_useTaper==1 ) + { + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, ux_np1 ); + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, ux_n ); + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, uy_np1 ); + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, uy_n ); + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, uz_np1 ); + TaperKernel::multiplyByTaperCoeff< EXEC_POLICY >( nodeManager.size(), taperCoeff, uz_n ); + } + + } void ElasticWaveEquationSEM::synchronizeUnknowns( real64 const & time_n, 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 8907f35cd76..b5f3ef3599d 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp @@ -157,6 +157,11 @@ class ElasticWaveEquationSEM : public WaveSolverBase void prepareNextTimestep( MeshLevel & mesh ); + /** + * @brief Get the minimum wavespeed on a mesh (S-wavespeed in the elastic case) + */ + virtual real32 getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /** * @brief Computes the minimum attenuation quality factor over all the mesh. This is useful for computing anelasticity coefficients, which * are usually global parameters diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/TaperKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/TaperKernel.hpp new file mode 100644 index 00000000000..c5b0f01c0a0 --- /dev/null +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/TaperKernel.hpp @@ -0,0 +1,158 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TaperKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_TAPERKERNEL_HPP_ +#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_TAPERKERNEL_HPP_ + +#include "WaveSolverUtils.hpp" + +namespace geos +{ + +using wsCoordType = WaveSolverUtils::wsCoordType; + +struct TaperKernel +{ + + /** + * @brief Compute coefficients for the taper layers + * @tparam EXEC_POLICY the execution policy + * @param[in] size the number of nodes + * @param[in] nodeCoords Coordinates of the nodes of the mesh (included interior degrees of freedom) + * @param[in] sizeT Taper thickness + * @param[in] dt time-step + * @param[in] vMin Min wavespeed (P-wavespeed for acoustic, S-wavespeed for elastic) + * @param[in] r desired reflectivity of the Taper + * @param[out] taperCoeff array which contains the taper coefficient on each node (which will be equal to 1 when we are outside of the + * taper layers) + */ + template< typename EXEC_POLICY > + static void + computeTaperCoeff( localIndex const size, + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + real32 const sizeT, + real32 const dt, + real32 const vMin, + real32 const r, + arrayView1d< real32 > const taperCoeff ) + { + + ///Seek the global maximum and minimum of the domain + RAJA::ReduceMin< parallelDeviceReduce, real32 > xMinGlobal( LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMin< parallelDeviceReduce, real32 > yMinGlobal( LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMin< parallelDeviceReduce, real32 > zMinGlobal( LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMax< parallelDeviceReduce, real32 > xMaxGlobal( -LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMax< parallelDeviceReduce, real32 > yMaxGlobal( -LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMax< parallelDeviceReduce, real32 > zMaxGlobal( -LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMin< parallelDeviceReduce, real32 > xMinInterior( LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMin< parallelDeviceReduce, real32 > yMinInterior( LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMin< parallelDeviceReduce, real32 > zMinInterior( LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMax< parallelDeviceReduce, real32 > xMaxInterior( -LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMax< parallelDeviceReduce, real32 > yMaxInterior( -LvArray::NumericLimits< real32 >::max ); + RAJA::ReduceMax< parallelDeviceReduce, real32 > zMaxInterior( -LvArray::NumericLimits< real32 >::max ); + + real32 xGlobalMin[3]{}; + real32 xGlobalMax[3]{}; + + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + xMinGlobal.min( nodeCoords[a][0] ); + yMinGlobal.min( nodeCoords[a][1] ); + zMinGlobal.min( nodeCoords[a][2] ); + xMaxGlobal.max( nodeCoords[a][0] ); + yMaxGlobal.max( nodeCoords[a][1] ); + zMaxGlobal.max( nodeCoords[a][2] ); + } ); + + xGlobalMin[0] = xMinGlobal.get(); + xGlobalMin[1] = yMinGlobal.get(); + xGlobalMin[2] = zMinGlobal.get(); + xGlobalMax[0] = xMaxGlobal.get(); + xGlobalMax[1] = yMaxGlobal.get(); + xGlobalMax[2] = zMaxGlobal.get(); + + + for( integer i=0; i<3; ++i ) + { + xGlobalMin[i] = MpiWrapper::min( xGlobalMin[i] ); + xGlobalMax[i] = MpiWrapper::max( xGlobalMax[i] ); + } + + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + real32 dist=0; + + real32 distXmin = LvArray::math::abs( nodeCoords[a][0]-xGlobalMin[0] ); + real32 distXmax = LvArray::math::abs( nodeCoords[a][0]-xGlobalMax[0] ); + real32 distYmin = LvArray::math::abs( nodeCoords[a][1]-xGlobalMin[1] ); + real32 distYmax = LvArray::math::abs( nodeCoords[a][1]-xGlobalMax[1] ); + real32 distZmax = LvArray::math::abs( nodeCoords[a][2]-xGlobalMax[2] ); + + dist=LvArray::math::min( distXmin, (LvArray::math::min( distXmax, LvArray::math::min( distYmin, (LvArray::math::min( distYmax, distZmax )))))); + + taperCoeff[a] = dist; + + + } ); + + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + + real32 dist = taperCoeff[a]; + + if( dist + static void + multiplyByTaperCoeff( localIndex const size, + arrayView1d< real32 const > const taperCoeff, + arrayView1d< real32 > const vector ) + { + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + vector[a] *= taperCoeff[a]; + } ); + + } + + +}; + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_TAPERERNEL_HPP_ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp index b55ca1f38d7..2114b824c12 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp @@ -207,6 +207,21 @@ WaveSolverBase::WaveSolverBase( const std::string & name, setSizedFromParent( 0 ). setDescription( "Element containing the receivers" ); + registerWrapper( viewKeyStruct::useTaperString(), &m_useTaper ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "Flag to apply taper" ); + + registerWrapper( viewKeyStruct::reflectivityCoeffString(), &m_reflectivityCoeff ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0.001 ). + setDescription( "Reflectivity coeff for taper" ); + + registerWrapper( viewKeyStruct::thicknessTaperString(), &m_thicknessTaper ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0.0 ). + setDescription( "Size for the taper layer " ); + registerWrapper( viewKeyStruct::slsReferenceAngularFrequenciesString(), &m_slsReferenceAngularFrequencies ). setInputFlag( InputFlags::OPTIONAL ). setSizedFromParent( 0 ). @@ -265,6 +280,10 @@ void WaveSolverBase::registerDataOnMesh( Group & meshBodies ) nodeCoords32[i][j] = X[i][j]; } } + + + nodeManager.registerField< fields::taperCoeff >( this->getName()); + } ); } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp index 2c319ca4f38..812fd73e2ed 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp @@ -119,6 +119,10 @@ class WaveSolverBase : public PhysicsSolverBase static constexpr char const * usePMLString() { return "usePML"; } static constexpr char const * parametersPMLString() { return "parametersPML"; } + static constexpr char const * useTaperString() {return "useTaper";} + static constexpr char const * reflectivityCoeffString() {return "reflectivityCoeff";} + static constexpr char const * thicknessTaperString() {return "thicknessTaper";} + static constexpr char const * receiverElemString() { return "receiverElem"; } static constexpr char const * receiverRegionString() { return "receiverRegion"; } static constexpr char const * freeSurfaceString() { return "FreeSurface"; } @@ -249,6 +253,11 @@ class WaveSolverBase : public PhysicsSolverBase DomainPartition & domain, bool const computeGradient ) = 0; + /** + * @brief Method to get the maximum wavespeed on a mesh (usually the P-wavespeed) + */ + virtual real32 getGlobalMinWavespeed( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) = 0; + virtual void registerDataOnMesh( Group & meshBodies ) override; @@ -320,6 +329,9 @@ class WaveSolverBase : public PhysicsSolverBase /// Flag to apply PML integer m_usePML; + ///Flag to use a taper + integer m_useTaper; + /// Flag to precompute the time-step /// usage: the time-step is computed then the code exit and you can /// copy paste the time-step inside the XML then deactivate the option @@ -370,6 +382,11 @@ class WaveSolverBase : public PhysicsSolverBase /// A set of target nodes IDs that will be handled by the current solver SortedArray< localIndex > m_solverTargetNodesSet; + /// Thickness of the Taper region, used to compute the damping profile + real32 m_thicknessTaper; + + real32 m_reflectivityCoeff; + /// Names of table functions for source wavelet (time dependency) array1d< string > m_sourceWaveletTableNames; @@ -399,6 +416,7 @@ class WaveSolverBase : public PhysicsSolverBase R1Tensor32 waveSpeedMaxXYZPML; }; + }; namespace fields @@ -411,6 +429,14 @@ DECLARE_FIELD( referencePosition32, NOPLOT, WRITE_AND_READ, "Copy of the referencePosition from NodeManager in 32 bits integer" ); +DECLARE_FIELD( taperCoeff, + "taperCoeff", + array1d< real32 >, + 1.0, + LEVEL_0, + WRITE_AND_READ, + "Array continaing the coefficients for the taper" ); + } } /* namespace geos */ diff --git a/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst b/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst new file mode 100644 index 00000000000..2e5e71e8b79 --- /dev/null +++ b/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst @@ -0,0 +1,43 @@ + + +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +attenuationType geos_WaveSolverUtils_AttenuationType none Flag to indicate which attenuation model to use: "none" for no attenuation, "sls\ for the standard-linear-solid (SLS) model (Fichtner, 2014). +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +dtSeismoTrace real64 0 Time step for output pressure at receivers +enableLifo integer 0 Set to 1 to enable LIFO storage feature +forward integer 1 Set to 1 to compute forward propagation +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +lifoOnDevice integer -80 Set the capacity of the lifo device storage (if negative, opposite of percentage of remaining memory) +lifoOnHost integer -80 Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory) +lifoSize integer 2147483647 Set the capacity of the lifo storage (should be the total number of buffers to store in the LIFO) +linearDASGeometry real64_array2d {{0}} Geometry parameters for a linear DAS fiber (dip, azimuth, gauge length) +linearDASSamples integer 5 Number of sample points to be used for strain integration when integrating the strain for the DAS signal +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise +receiverCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the receivers +reflectivityCoeff real32 0.001 Reflectivity coeff for taper +rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default +saveFields integer 0 Set to 1 to save fields during forward and restore them during backward +shotIndex integer 0 Set the current shot for temporary files +slsAnelasticityCoefficients real32_array {0} Anelasticity coefficients for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding reference frequencies must be provided. +slsReferenceAngularFrequencies real32_array {0} Reference angular frequencies (omega) for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding anelasticity coefficients must be provided. +sourceCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the sources +targetRegions groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +thicknessMaxXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +thicknessMinXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +timeSourceDelay real32 -1 Source time delay (1 / f0 by default) +timeSourceFrequency real32 0 Central frequency for the time source +useDAS geos_WaveSolverUtils_DASType none Flag to indicate if DAS data will be modeled, and which DAS type to use: "none" to deactivate DAS, "strainIntegration" for strain integration, "dipole" for displacement difference +useTaper integer 0 Flag to apply taper +writeLinearSystem integer 0 Write matrix, rhs, solution to screen ( = 1) or file ( = 2). +xMaxTaper R1Tensor32 {0,0,0} Maximal coordinates for taper (right,top,back) +xMinTaper R1Tensor32 {0,0,0} Minimal coordinates for taper (left,bottom,front) +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/schema/docs/AcousticSEM.rst b/src/coreComponents/schema/docs/AcousticSEM.rst new file mode 100644 index 00000000000..2e5e71e8b79 --- /dev/null +++ b/src/coreComponents/schema/docs/AcousticSEM.rst @@ -0,0 +1,43 @@ + + +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +attenuationType geos_WaveSolverUtils_AttenuationType none Flag to indicate which attenuation model to use: "none" for no attenuation, "sls\ for the standard-linear-solid (SLS) model (Fichtner, 2014). +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +dtSeismoTrace real64 0 Time step for output pressure at receivers +enableLifo integer 0 Set to 1 to enable LIFO storage feature +forward integer 1 Set to 1 to compute forward propagation +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +lifoOnDevice integer -80 Set the capacity of the lifo device storage (if negative, opposite of percentage of remaining memory) +lifoOnHost integer -80 Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory) +lifoSize integer 2147483647 Set the capacity of the lifo storage (should be the total number of buffers to store in the LIFO) +linearDASGeometry real64_array2d {{0}} Geometry parameters for a linear DAS fiber (dip, azimuth, gauge length) +linearDASSamples integer 5 Number of sample points to be used for strain integration when integrating the strain for the DAS signal +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise +receiverCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the receivers +reflectivityCoeff real32 0.001 Reflectivity coeff for taper +rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default +saveFields integer 0 Set to 1 to save fields during forward and restore them during backward +shotIndex integer 0 Set the current shot for temporary files +slsAnelasticityCoefficients real32_array {0} Anelasticity coefficients for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding reference frequencies must be provided. +slsReferenceAngularFrequencies real32_array {0} Reference angular frequencies (omega) for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding anelasticity coefficients must be provided. +sourceCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the sources +targetRegions groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +thicknessMaxXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +thicknessMinXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +timeSourceDelay real32 -1 Source time delay (1 / f0 by default) +timeSourceFrequency real32 0 Central frequency for the time source +useDAS geos_WaveSolverUtils_DASType none Flag to indicate if DAS data will be modeled, and which DAS type to use: "none" to deactivate DAS, "strainIntegration" for strain integration, "dipole" for displacement difference +useTaper integer 0 Flag to apply taper +writeLinearSystem integer 0 Write matrix, rhs, solution to screen ( = 1) or file ( = 2). +xMaxTaper R1Tensor32 {0,0,0} Maximal coordinates for taper (right,top,back) +xMinTaper R1Tensor32 {0,0,0} Minimal coordinates for taper (left,bottom,front) +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/schema/docs/AcousticVTISEM.rst b/src/coreComponents/schema/docs/AcousticVTISEM.rst new file mode 100644 index 00000000000..2e5e71e8b79 --- /dev/null +++ b/src/coreComponents/schema/docs/AcousticVTISEM.rst @@ -0,0 +1,43 @@ + + +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +attenuationType geos_WaveSolverUtils_AttenuationType none Flag to indicate which attenuation model to use: "none" for no attenuation, "sls\ for the standard-linear-solid (SLS) model (Fichtner, 2014). +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +dtSeismoTrace real64 0 Time step for output pressure at receivers +enableLifo integer 0 Set to 1 to enable LIFO storage feature +forward integer 1 Set to 1 to compute forward propagation +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +lifoOnDevice integer -80 Set the capacity of the lifo device storage (if negative, opposite of percentage of remaining memory) +lifoOnHost integer -80 Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory) +lifoSize integer 2147483647 Set the capacity of the lifo storage (should be the total number of buffers to store in the LIFO) +linearDASGeometry real64_array2d {{0}} Geometry parameters for a linear DAS fiber (dip, azimuth, gauge length) +linearDASSamples integer 5 Number of sample points to be used for strain integration when integrating the strain for the DAS signal +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise +receiverCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the receivers +reflectivityCoeff real32 0.001 Reflectivity coeff for taper +rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default +saveFields integer 0 Set to 1 to save fields during forward and restore them during backward +shotIndex integer 0 Set the current shot for temporary files +slsAnelasticityCoefficients real32_array {0} Anelasticity coefficients for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding reference frequencies must be provided. +slsReferenceAngularFrequencies real32_array {0} Reference angular frequencies (omega) for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding anelasticity coefficients must be provided. +sourceCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the sources +targetRegions groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +thicknessMaxXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +thicknessMinXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +timeSourceDelay real32 -1 Source time delay (1 / f0 by default) +timeSourceFrequency real32 0 Central frequency for the time source +useDAS geos_WaveSolverUtils_DASType none Flag to indicate if DAS data will be modeled, and which DAS type to use: "none" to deactivate DAS, "strainIntegration" for strain integration, "dipole" for displacement difference +useTaper integer 0 Flag to apply taper +writeLinearSystem integer 0 Write matrix, rhs, solution to screen ( = 1) or file ( = 2). +xMaxTaper R1Tensor32 {0,0,0} Maximal coordinates for taper (right,top,back) +xMinTaper R1Tensor32 {0,0,0} Minimal coordinates for taper (left,bottom,front) +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst b/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst new file mode 100644 index 00000000000..2e5e71e8b79 --- /dev/null +++ b/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst @@ -0,0 +1,43 @@ + + +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== +attenuationType geos_WaveSolverUtils_AttenuationType none Flag to indicate which attenuation model to use: "none" for no attenuation, "sls\ for the standard-linear-solid (SLS) model (Fichtner, 2014). +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +dtSeismoTrace real64 0 Time step for output pressure at receivers +enableLifo integer 0 Set to 1 to enable LIFO storage feature +forward integer 1 Set to 1 to compute forward propagation +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +lifoOnDevice integer -80 Set the capacity of the lifo device storage (if negative, opposite of percentage of remaining memory) +lifoOnHost integer -80 Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory) +lifoSize integer 2147483647 Set the capacity of the lifo storage (should be the total number of buffers to store in the LIFO) +linearDASGeometry real64_array2d {{0}} Geometry parameters for a linear DAS fiber (dip, azimuth, gauge length) +linearDASSamples integer 5 Number of sample points to be used for strain integration when integrating the strain for the DAS signal +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise +receiverCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the receivers +reflectivityCoeff real32 0.001 Reflectivity coeff for taper +rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default +saveFields integer 0 Set to 1 to save fields during forward and restore them during backward +shotIndex integer 0 Set the current shot for temporary files +slsAnelasticityCoefficients real32_array {0} Anelasticity coefficients for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding reference frequencies must be provided. +slsReferenceAngularFrequencies real32_array {0} Reference angular frequencies (omega) for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding anelasticity coefficients must be provided. +sourceCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the sources +targetRegions groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +thicknessMaxXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +thicknessMinXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +timeSourceDelay real32 -1 Source time delay (1 / f0 by default) +timeSourceFrequency real32 0 Central frequency for the time source +useDAS geos_WaveSolverUtils_DASType none Flag to indicate if DAS data will be modeled, and which DAS type to use: "none" to deactivate DAS, "strainIntegration" for strain integration, "dipole" for displacement difference +useTaper integer 0 Flag to apply taper +writeLinearSystem integer 0 Write matrix, rhs, solution to screen ( = 1) or file ( = 2). +xMaxTaper R1Tensor32 {0,0,0} Maximal coordinates for taper (right,top,back) +xMinTaper R1Tensor32 {0,0,0} Minimal coordinates for taper (left,bottom,front) +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +============================== ==================================== ========== ======================================================================================================================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/schema/docs/ElasticSEM.rst b/src/coreComponents/schema/docs/ElasticSEM.rst new file mode 100644 index 00000000000..1f3558932eb --- /dev/null +++ b/src/coreComponents/schema/docs/ElasticSEM.rst @@ -0,0 +1,46 @@ + + +============================== ==================================== ============= ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +============================== ==================================== ============= ======================================================================================================================================================================================================================================================================================================================== +attenuationType geos_WaveSolverUtils_AttenuationType none Flag to indicate which attenuation model to use: "none" for no attenuation, "sls\ for the standard-linear-solid (SLS) model (Fichtner, 2014). +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +dtSeismoTrace real64 0 Time step for output pressure at receivers +enableLifo integer 0 Set to 1 to enable LIFO storage feature +forward integer 1 Set to 1 to compute forward propagation +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +lifoOnDevice integer -80 Set the capacity of the lifo device storage (if negative, opposite of percentage of remaining memory) +lifoOnHost integer -80 Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory) +lifoSize integer 2147483647 Set the capacity of the lifo storage (should be the total number of buffers to store in the LIFO) +linearDASGeometry real64_array2d {{0}} Geometry parameters for a linear DAS fiber (dip, azimuth, gauge length) +linearDASSamples integer 5 Number of sample points to be used for strain integration when integrating the strain for the DAS signal +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise +receiverCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the receivers +reflectivityCoeff real32 0.001 Reflectivity coeff for taper +rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default +saveFields integer 0 Set to 1 to save fields during forward and restore them during backward +shotIndex integer 0 Set the current shot for temporary files +slsAnelasticityCoefficients real32_array {0} Anelasticity coefficients for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding reference frequencies must be provided. +slsReferenceAngularFrequencies real32_array {0} Reference angular frequencies (omega) for the standard-linear-solid (SLS) anelasticity.The default value is { }, corresponding to no attenuation. An array with the corresponding anelasticity coefficients must be provided. +sourceCoordinates real64_array2d {{0}} Coordinates (x,y,z) of the sources +sourceForce R1Tensor {0,0,0} Force of the source: 3 real values for a vector source, and 6 real values for a tensor source (in Voigt notation).The default value is { 0, 0, 0 } (no net force). +sourceMoment R2SymTensor {1,1,1,0,0,0} Moment of the source: 6 real values describing a symmetric tensor in Voigt notation.The default value is { 1, 1, 1, 0, 0, 0 } (diagonal moment, corresponding to a pure explosion). +targetRegions groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +thicknessMaxXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +thicknessMinXYZTaper R1Tensor32 {0,0,0} Size for the taper layer (left,bottom,front) +timeSourceDelay real32 -1 Source time delay (1 / f0 by default) +timeSourceFrequency real32 0 Central frequency for the time source +useDAS geos_WaveSolverUtils_DASType none Flag to indicate if DAS data will be modeled, and which DAS type to use: "none" to deactivate DAS, "strainIntegration" for strain integration, "dipole" for displacement difference +useTaper integer 0 Flag to apply taper +useVTI integer 0 Flag to apply VTI anisotropy. The default is to use isotropic physic. +writeLinearSystem integer 0 Write matrix, rhs, solution to screen ( = 1) or file ( = 2). +xMaxTaper R1Tensor32 {0,0,0} Maximal coordinates for taper (right,top,back) +xMinTaper R1Tensor32 {0,0,0} Minimal coordinates for taper (left,bottom,front) +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +============================== ==================================== ============= ======================================================================================================================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index e130b1e56d4..0a62161b03c 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -2404,6 +2404,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2420,6 +2422,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2428,6 +2432,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2489,6 +2495,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2505,6 +2513,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2513,6 +2523,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2564,6 +2576,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2580,6 +2594,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2588,6 +2604,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -3029,6 +3047,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -3045,6 +3065,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -3053,6 +3075,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -3104,6 +3128,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -3124,6 +3150,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -3132,6 +3160,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + diff --git a/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt b/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt index c1296f06b10..75412c2a24d 100644 --- a/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt @@ -1,6 +1,7 @@ # Specify list of tests set( gtest_geosx_tests testWavePropagation.cpp + testWavePropagationTaper.cpp testWavePropagationQ2.cpp testWavePropagationElasticFirstOrder.cpp testWavePropagationDAS.cpp diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationTaper.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationTaper.cpp new file mode 100644 index 00000000000..f5721e4c19e --- /dev/null +++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationTaper.cpp @@ -0,0 +1,244 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 Total, S.A + * Copyright (c) 2020- GEOSX Contributors + * All right 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 pReceivers = propagator->getReference< array2d< real32 > >( AcousticWaveEquationSEM::viewKeyStruct::pressureNp1AtReceiversString() ).toView(); + + // move it to CPU, if needed + pReceivers.move( hostMemorySpace, false ); + + // check number of seismos and trace length + ASSERT_EQ( pReceivers.size( 1 ), 10 ); + ASSERT_EQ( pReceivers.size( 0 ), 11 ); + + // check seismo content. The pressure 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( pReceivers[i][8] ) <= 0.0000001 ); + } + } + // check again the seismo content. + for( int i = 0; i < 11; i++ ) + { + if( i > 0 ) + { + ASSERT_TRUE( std::abs( pReceivers[i][8] ) <= 0.0000001 ); + } + } +} + +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; +}