From ece2ae97086d1ccc98da37d71c3f7c7340d2222b Mon Sep 17 00:00:00 2001 From: tbeltzun <129868353+tbeltzun@users.noreply.github.com> Date: Thu, 14 Dec 2023 04:11:50 +0100 Subject: [PATCH] fix elastic seismos bug (#2849) The rewrite of #2729 introduced a bug in the elastic solver, where the derived class had members shadowing the parent class variables. The cleanup was done in #2557 but was somehow missed when cherry-picking commits to create #2729. This PR: - fixes elastic solver bug (member variables shadowing), reported by @shenchengyi, reproducer using `inputFiles/wavePropagation/elas3D_DAS_smoke.xml`: empty seismos time and sampled value; - make `initTrace` consistent with `writeSeismoTrace` i.e. don't touch the file if `outputSeismoTrace` is unset; - remove redundant / unused `ElasticWaveEquationSEM::getNumNodesPerElem`; - code style consistency across `wavePropagation` solvers files (e.g. use `e` for element); - use `pow(..., n)` patterns for clarity. --- integratedTests | 2 +- .../AcousticFirstOrderWaveEquationSEM.cpp | 4 +- ...cousticFirstOrderWaveEquationSEMKernel.hpp | 15 +++-- .../AcousticVTIWaveEquationSEM.cpp | 4 +- .../AcousticVTIWaveEquationSEMKernel.hpp | 9 ++- .../AcousticWaveEquationSEM.cpp | 10 ++- .../AcousticWaveEquationSEMKernel.hpp | 16 ++--- .../ElasticFirstOrderWaveEquationSEM.cpp | 2 - ...ElasticFirstOrderWaveEquationSEMKernel.hpp | 10 +-- .../ElasticWaveEquationSEM.cpp | 61 +++---------------- .../ElasticWaveEquationSEM.hpp | 17 ------ .../ElasticWaveEquationSEMKernel.hpp | 27 ++++---- .../wavePropagation/WaveSolverUtils.hpp | 8 ++- 13 files changed, 57 insertions(+), 128 deletions(-) diff --git a/integratedTests b/integratedTests index b08da69e4e0..5ffebf05d09 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit b08da69e4e0f57b744b4f5496a829322b8fc3820 +Subproject commit 5ffebf05d093b7492328e782c68c1612e6ce14af diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp index c7ec3d9d330..fc8fb56b5a2 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp @@ -210,7 +210,6 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev { using FE_TYPE = TYPEOFREF( finiteElement ); - constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); acousticFirstOrderWaveEquationSEMKernels:: @@ -218,7 +217,6 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev launch< EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), regionIndex, - numNodesPerElem, numFacesPerElem, X, elemGhostRank, @@ -343,7 +341,7 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro } ); } ); - WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp index ac793524a2a..9347166af08 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp @@ -39,7 +39,6 @@ struct PrecomputeSourceAndReceiverKernel * @tparam EXEC_POLICY execution policy * @tparam FE_TYPE finite element type * @param[in] size the number of cells in the subRegion - * @param[in] numNodesPerElem number of nodes per element * @param[in] numFacesPerElem number of faces per element * @param[in] nodeCoords coordinates of the nodes * @param[in] elemGhostRank rank of the ghost element @@ -68,7 +67,6 @@ struct PrecomputeSourceAndReceiverKernel static void launch( localIndex const size, localIndex const regionIndex, - localIndex const numNodesPerElem, localIndex const numFacesPerElem, arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, arrayView1d< integer const > const elemGhostRank, @@ -95,6 +93,7 @@ struct PrecomputeSourceAndReceiverKernel real32 const timeSourceDelay, localIndex const rickerOrder ) { + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { @@ -133,7 +132,7 @@ struct PrecomputeSourceAndReceiverKernel sourceIsAccessible[isrc] = 1; sourceElem[isrc] = k; sourceRegion[isrc] = regionIndex; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -181,7 +180,7 @@ struct PrecomputeSourceAndReceiverKernel rcvElem[ircv] = k; receiverRegion[ircv] = regionIndex; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -227,26 +226,26 @@ struct MassMatrixKernel arrayView1d< real32 > const mass ) { - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) { constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; - real32 const invC2 = 1.0 / ( density[k] * velocity[k] * velocity[k] ); + real32 const invC2 = 1.0 / ( density[e] * pow( velocity[e], 2 ) ); real64 xLocal[ numNodesPerElem ][ 3 ]; for( localIndex a = 0; a < numNodesPerElem; ++a ) { for( localIndex i = 0; i < 3; ++i ) { - xLocal[a][i] = nodeCoords( elemsToNodes( k, a ), i ); + xLocal[a][i] = nodeCoords( elemsToNodes( e, a ), i ); } } for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q ) { real32 const localIncrement = invC2 * m_finiteElement.computeMassTerm( q, xLocal ); - RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement ); + RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( e, q )], localIncrement ); } } ); // end loop over element } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp index 8922a92e8cd..28851545659 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp @@ -168,14 +168,12 @@ void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & me { using FE_TYPE = TYPEOFREF( finiteElement ); - constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); acousticVTIWaveEquationSEMKernels:: PrecomputeSourceAndReceiverKernel:: launch< EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), - numNodesPerElem, numFacesPerElem, X32, elemGhostRank, @@ -314,7 +312,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); - WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); } void AcousticVTIWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain ) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp index bd062c38fa2..52434125c57 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp @@ -39,7 +39,6 @@ struct PrecomputeSourceAndReceiverKernel * @tparam EXEC_POLICY execution policy * @tparam FE_TYPE finite element type * @param[in] size the number of cells in the subRegion - * @param[in] numNodesPerElem number of nodes per element * @param[in] numFacesPerElem number of faces per element * @param[in] nodeCoords coordinates of the nodes * @param[in] elemGhostRank the ghost ranks @@ -64,7 +63,6 @@ struct PrecomputeSourceAndReceiverKernel template< typename EXEC_POLICY, typename FE_TYPE > static void launch( localIndex const size, - localIndex const numNodesPerElem, localIndex const numFacesPerElem, arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, arrayView1d< integer const > const elemGhostRank, @@ -87,6 +85,7 @@ struct PrecomputeSourceAndReceiverKernel real32 const timeSourceDelay, localIndex const rickerOrder ) { + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { @@ -123,7 +122,7 @@ struct PrecomputeSourceAndReceiverKernel coordsOnRefElem ); sourceIsAccessible[isrc] = 1; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -170,7 +169,7 @@ struct PrecomputeSourceAndReceiverKernel receiverIsLocal[ircv] = 1; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -220,7 +219,7 @@ struct MassMatrixKernel constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; - real32 const invC2 = 1.0 / ( velocity[e] * velocity[e] ); + real32 const invC2 = 1.0 / pow( velocity[e], 2 ); real64 xLocal[ numNodesPerElem ][ 3 ]; for( localIndex a = 0; a < numNodesPerElem; ++a ) { diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp index 6cbadc9ec14..eb4ed4d8329 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp @@ -172,7 +172,6 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, { using FE_TYPE = TYPEOFREF( finiteElement ); - constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); { @@ -181,7 +180,6 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, PrecomputeSourceAndReceiverKernel:: launch< EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), - numNodesPerElem, numFacesPerElem, X32, elemGhostRank, @@ -323,7 +321,7 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); - WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); } @@ -853,7 +851,7 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n, EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() ); - int const maxCycle = int(round( maxTime/dt )); + int const maxCycle = int(round( maxTime / dt )); if( computeGradient && cycleNumber < maxCycle ) { @@ -962,8 +960,8 @@ real64 AcousticWaveEquationSEM::explicitStepInternal( real64 const & time_n, EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() ); - integer const cycleForSource = int(round( -minTime/dt + cycleNumber )); - //std::cout<<"cycle GEOSX = "< static void launch( localIndex const size, - localIndex const numNodesPerElem, localIndex const numFacesPerElem, arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, arrayView1d< integer const > const elemGhostRank, @@ -81,9 +79,11 @@ struct PrecomputeSourceAndReceiverKernel real32 const timeSourceDelay, localIndex const rickerOrder ) { + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { + real64 const center[3] = { elemCenter[k][0], elemCenter[k][1], elemCenter[k][2] }; @@ -117,7 +117,7 @@ struct PrecomputeSourceAndReceiverKernel coordsOnRefElem ); sourceIsAccessible[isrc] = 1; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -164,7 +164,7 @@ struct PrecomputeSourceAndReceiverKernel receiverIsLocal[ircv] = 1; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -216,7 +216,7 @@ struct MassMatrixKernel constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; - real32 const invC2 = 1.0 / ( density[e] * velocity[e] * velocity[e] ); + real32 const invC2 = 1.0 / ( density[e] * pow( velocity[e], 2 ) ); real64 xLocal[ numNodesPerElem ][ 3 ]; for( localIndex a = 0; a < numNodesPerElem; ++a ) { @@ -803,9 +803,9 @@ class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, { m_finiteElementSpace.template computeStiffnessTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { - real32 invDensity = 1./m_density[k]; - real32 const localIncrement = invDensity*val*m_p_n[m_elemsToNodes[k][j]]; - RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector[m_elemsToNodes[k][i]], localIncrement ); + real32 invDensity = 1.0 / m_density[k]; + real32 const localIncrement = invDensity*val*m_p_n[m_elemsToNodes( k, j )]; + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector[m_elemsToNodes( k, i )], localIncrement ); } ); } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp index 1e4364f840c..29c61faa2e5 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp @@ -260,7 +260,6 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve { using FE_TYPE = TYPEOFREF( finiteElement ); - constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); elasticFirstOrderWaveEquationSEMKernels:: @@ -268,7 +267,6 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve launch< EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), regionIndex, - numNodesPerElem, numFacesPerElem, X, elemGhostRank, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp index c4414f17c0e..262fe06824e 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp @@ -62,7 +62,6 @@ struct PrecomputeSourceAndReceiverKernel static void launch( localIndex const size, localIndex const regionIndex, - localIndex const numNodesPerElem, localIndex const numFacesPerElem, arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, arrayView1d< integer const > const elemGhostRank, @@ -89,6 +88,7 @@ struct PrecomputeSourceAndReceiverKernel real32 const timeSourceDelay, localIndex const rickerOrder ) { + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { @@ -126,7 +126,7 @@ struct PrecomputeSourceAndReceiverKernel sourceIsAccessible[isrc] = 1; sourceElem[isrc] = k; sourceRegion[isrc] = regionIndex; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -174,7 +174,7 @@ struct PrecomputeSourceAndReceiverKernel rcvElem[ircv] = k; receiverRegion[ircv] = regionIndex; - real64 Ntest[FE_TYPE::numNodes]; + real64 Ntest[numNodesPerElem]; FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -377,8 +377,8 @@ struct StressComputation } } - mu[k] = density[k] * velocityVs[k] * velocityVs[k]; - lambda[k] = density[k] * velocityVp[k] * velocityVp[k] - 2.0*mu[k]; + mu[k] = density[k] * pow( velocityVs[k], 2 ); + lambda[k] = density[k] * pow( velocityVp[k], 2 ) - 2.0*mu[k]; real32 uelemxx[numNodesPerElem] = {0.0}; real32 uelemyy[numNodesPerElem] = {0.0}; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp index 2de2021b306..b52daa13cb1 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp @@ -88,49 +88,6 @@ ElasticWaveEquationSEM::~ElasticWaveEquationSEM() // TODO Auto-generated destructor stub } -localIndex ElasticWaveEquationSEM::getNumNodesPerElem() -{ - DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); - - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - - FiniteElementDiscretizationManager const & - feDiscretizationManager = numericalMethodManager.getFiniteElementDiscretizationManager(); - - FiniteElementDiscretization const * const - feDiscretization = feDiscretizationManager.getGroupPointer< FiniteElementDiscretization >( m_discretizationName ); - GEOS_THROW_IF( feDiscretization == nullptr, - getDataContext() << ": FE discretization not found: " << m_discretizationName, - InputError ); - - localIndex numNodesPerElem = 0; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), - [&]( string const &, - MeshLevel const & mesh, - arrayView1d< string const > const & regionNames ) - { - ElementRegionManager const & elemManager = mesh.getElemManager(); - elemManager.forElementRegions( regionNames, - [&] ( localIndex const, - ElementRegionBase const & elemRegion ) - { - elemRegion.forElementSubRegions( [&]( ElementSubRegionBase const & elementSubRegion ) - { - finiteElement::FiniteElementBase const & - fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); - localIndex const numSupportPoints = fe.getNumSupportPoints(); - if( numSupportPoints > numNodesPerElem ) - { - numNodesPerElem = numSupportPoints; - } - } ); - } ); - - - } ); - return numNodesPerElem; -} - void ElasticWaveEquationSEM::initializePreSubGroups() { @@ -234,7 +191,7 @@ void ElasticWaveEquationSEM::postProcessInput() { m_nsamplesSeismoTrace = 0; } - localIndex const nsamples = int(maxTime/dt) + 1; + localIndex const nsamples = int(maxTime / dt) + 1; localIndex const numSourcesGlobal = m_sourceCoordinates.size( 0 ); m_sourceIsAccessible.resize( numSourcesGlobal ); @@ -353,9 +310,9 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, } ); } -void ElasticWaveEquationSEM::computeDAS ( arrayView2d< real32 > const xCompRcv, - arrayView2d< real32 > const yCompRcv, - arrayView2d< real32 > const zCompRcv ) +void ElasticWaveEquationSEM::computeDAS( arrayView2d< real32 > const xCompRcv, + arrayView2d< real32 > const yCompRcv, + arrayView2d< real32 > const zCompRcv ) { arrayView2d< real64 const > const linearDASGeometry = m_linearDASGeometry.toViewConst(); @@ -399,9 +356,9 @@ void ElasticWaveEquationSEM::computeDAS ( arrayView2d< real32 > const xCompRcv, { // store strain data in the z-component of the receiver (copied to x after resize) zCompRcv[iSample][ircv] = - cd * ca * ( xCompRcv[iSample][numReceiversGlobal+ircv] - xCompRcv[iSample][ircv] ) - + cd * sa * ( yCompRcv[iSample][numReceiversGlobal+ircv] - yCompRcv[iSample][ircv] ) - + sd * ( zCompRcv[iSample][numReceiversGlobal+ircv] - zCompRcv[iSample][ircv] ); + cd * ca * ( xCompRcv[iSample][numReceiversGlobal+ircv] - xCompRcv[iSample][ircv] ) + + cd * sa * ( yCompRcv[iSample][numReceiversGlobal+ircv] - yCompRcv[iSample][ircv] ) + + sd * ( zCompRcv[iSample][numReceiversGlobal+ircv] - zCompRcv[iSample][ircv] ); zCompRcv[iSample][ircv] /= linearDASGeometry[ircv][2]; } @@ -544,8 +501,8 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); - WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); - WaveSolverUtils::initTrace( "dasTraceReceiver", getName(), m_linearDASGeometry.size( 0 ), m_receiverIsLocal ); + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); + WaveSolverUtils::initTrace( "dasTraceReceiver", getName(), m_outputSeismoTrace, m_linearDASGeometry.size( 0 ), m_receiverIsLocal ); } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp index 119b40da706..0cb7afd3d6f 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp @@ -171,11 +171,6 @@ class ElasticWaveEquationSEM : public WaveSolverBase */ virtual void applyPML( real64 const time, DomainPartition & domain ) override; - localIndex getNumNodesPerElem(); - - /// Indices of the nodes (in the right order) for each source point - array2d< localIndex > m_sourceNodeIds; - /// Constant part of the source for the nodes listed in m_sourceNodeIds in x-direction array2d< real64 > m_sourceConstantsx; @@ -185,18 +180,6 @@ class ElasticWaveEquationSEM : public WaveSolverBase /// Constant part of the source for the nodes listed in m_sourceNodeIds in x-direction array2d< real64 > m_sourceConstantsz; - /// Flag that indicates whether the source is accessible or not to the MPI rank - array1d< localIndex > m_sourceIsAccessible; - - /// Indices of the element nodes (in the right order) for each receiver point - array2d< localIndex > m_receiverNodeIds; - - /// Basis function evaluated at the receiver for the nodes listed in m_receiverNodeIds - array2d< real64 > m_receiverConstants; - - /// Flag that indicates whether the receiver is local or not to the MPI rank - array1d< localIndex > m_receiverIsLocal; - /// Displacement_np1 at the receiver location for each time step for each receiver (x-component) array2d< real32 > m_displacementXNp1AtReceivers; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp index ebdd56542f4..609f97bd245 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp @@ -90,12 +90,11 @@ struct PrecomputeSourceAndReceiverKernel R1Tensor const sourceForce, R2SymTensor const sourceMoment ) { + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { - constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; - real64 const center[3] = { elemCenter[k][0], elemCenter[k][1], elemCenter[k][2] }; @@ -141,8 +140,8 @@ struct PrecomputeSourceAndReceiverKernel coordsOnRefElem ); sourceIsAccessible[isrc] = 1; - real64 N[FE_TYPE::numNodes]; - real64 gradN[FE_TYPE::numNodes][3]; + real64 N[numNodesPerElem]; + real64 gradN[numNodesPerElem][3]; FE_TYPE::calcN( coordsOnRefElem, N ); FE_TYPE::calcGradN( coordsOnRefElem, xLocal, gradN ); R2SymTensor moment = sourceMoment; @@ -169,7 +168,6 @@ struct PrecomputeSourceAndReceiverKernel } } // end loop over all sources - // Step 2: locate the receivers, and precompute the receiver term /// loop over all the receivers that haven't been found yet @@ -196,11 +194,8 @@ struct PrecomputeSourceAndReceiverKernel elemsToNodes[k], nodeCoords, coordsOnRefElem ); - receiverIsLocal[ircv] = 1; - real64 Ntest[numNodesPerElem]; - FE_TYPE::calcN( coordsOnRefElem, Ntest ); for( localIndex a = 0; a < numNodesPerElem; ++a ) @@ -478,8 +473,8 @@ class ExplicitElasticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, stack.xLocal[ a ][ i ] = m_nodeCoords[ nodeIndex ][ i ]; } } - stack.mu = m_density[k] * m_velocityVs[k] * m_velocityVs[k]; - stack.lambda = m_density[k] *m_velocityVp[k] * m_velocityVp[k] - 2.0*stack.mu; + stack.mu = m_density[k] * pow( m_velocityVs[k], 2 ); + stack.lambda = m_density[k] * pow( m_velocityVp[k], 2 ) - 2.0 * stack.mu; } /** @@ -508,13 +503,13 @@ class ExplicitElasticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, real32 const Ryz_ij = val*(stack.lambda*J[p][1]*J[r][2]+stack.mu*J[p][2]*J[r][1]); real32 const Rzy_ij = val*(stack.mu*J[p][1]*J[r][2]+stack.lambda*J[p][2]*J[r][1]); - real32 const localIncrementx = (Rxx_ij * m_ux_n[m_elemsToNodes[k][j]] + Rxy_ij*m_uy_n[m_elemsToNodes[k][j]] + Rxz_ij*m_uz_n[m_elemsToNodes[k][j]]); - real32 const localIncrementy = (Ryx_ij * m_ux_n[m_elemsToNodes[k][j]] + Ryy_ij*m_uy_n[m_elemsToNodes[k][j]] + Ryz_ij*m_uz_n[m_elemsToNodes[k][j]]); - real32 const localIncrementz = (Rzx_ij * m_ux_n[m_elemsToNodes[k][j]] + Rzy_ij*m_uy_n[m_elemsToNodes[k][j]] + Rzz_ij*m_uz_n[m_elemsToNodes[k][j]]); + real32 const localIncrementx = (Rxx_ij * m_ux_n[m_elemsToNodes( k, j )] + Rxy_ij*m_uy_n[m_elemsToNodes( k, j )] + Rxz_ij*m_uz_n[m_elemsToNodes( k, j )]); + real32 const localIncrementy = (Ryx_ij * m_ux_n[m_elemsToNodes( k, j )] + Ryy_ij*m_uy_n[m_elemsToNodes( k, j )] + Ryz_ij*m_uz_n[m_elemsToNodes( k, j )]); + real32 const localIncrementz = (Rzx_ij * m_ux_n[m_elemsToNodes( k, j )] + Rzy_ij*m_uy_n[m_elemsToNodes( k, j )] + Rzz_ij*m_uz_n[m_elemsToNodes( k, j )]); - RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorx[m_elemsToNodes[k][i]], localIncrementx ); - RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectory[m_elemsToNodes[k][i]], localIncrementy ); - RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorz[m_elemsToNodes[k][i]], localIncrementz ); + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorx[m_elemsToNodes( k, i )], localIncrementx ); + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectory[m_elemsToNodes( k, i )], localIncrementy ); + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorz[m_elemsToNodes( k, i )], localIncrementz ); } ); } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp index ca52e3af4c0..1e856597658 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp @@ -80,9 +80,12 @@ struct WaveSolverUtils */ static void initTrace( char const * prefix, string const & name, + bool const outputSeismoTrace, localIndex const nReceivers, arrayView1d< localIndex const > const receiverIsLocal ) { + if( !outputSeismoTrace ) return; + string const outputDir = OutputBase::getOutputDirectory(); RAJA::ReduceSum< ReducePolicy< serialPolicy >, localIndex > count( 0 ); @@ -101,7 +104,7 @@ struct WaveSolverUtils } /** - * @brief Convenient helper for 3D vectors calling 3 times the scalar version. + * @brief Convenient helper for 3D vectors calling 3 times the scalar version with only the sampled variable argument changed. */ static void writeSeismoTraceVector( char const * prefix, string const & name, @@ -140,6 +143,7 @@ struct WaveSolverUtils std::ofstream f( fn, std::ios::app ); if( f ) { + GEOS_LOG_RANK( GEOS_FMT( "Append to seismo trace file {}", fn ) ); for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample ) { // index - time - value @@ -208,7 +212,7 @@ struct WaveSolverUtils arrayView2d< real32 const > const var_n, arrayView2d< real32 > varAtReceivers ) { - real64 const time_np1 = time_n+dt; + real64 const time_np1 = time_n + dt; real32 const a1 = dt < epsilonLoc ? 1.0 : (time_np1 - timeSeismo) / dt; real32 const a2 = 1.0 - a1;