From e706287a58f94c671b8485dd0c2453c436d91719 Mon Sep 17 00:00:00 2001 From: Bertrand Thierry Date: Tue, 10 Dec 2024 15:33:33 -0600 Subject: [PATCH] conflit solving + update solvers --- ...TIFletcherAdjointWaveEquationSEMKernel.hpp | 57 +-- .../AcousticVTIFletcherWaveEquationSEM.cpp | 349 +++++++++++------ .../AcousticVTIFletcherWaveEquationSEM.hpp | 24 +- ...ousticVTIFletcherWaveEquationSEMKernel.hpp | 45 ++- ...icVTIZhangAdjointWaveEquationSEMKernel.hpp | 49 ++- .../AcousticVTIZhangWaveEquationSEM.cpp | 361 ++++++++++-------- .../AcousticVTIZhangWaveEquationSEM.hpp | 33 +- .../AcousticVTIZhangWaveEquationSEMKernel.hpp | 46 +-- 8 files changed, 535 insertions(+), 429 deletions(-) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherAdjointWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherAdjointWaveEquationSEMKernel.hpp index d667d087b65..beeffaa1363 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherAdjointWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherAdjointWaveEquationSEMKernel.hpp @@ -14,7 +14,7 @@ */ /** - * @file AcousticVTIFletcherAdjointWaveEquationSEMKernel.hpp + * @file AcousticVTIFletcherWaveEquationSEMKernel.hpp */ #ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICVTIFLETCHERADJOINTWAVEEQUATIONSEMKERNEL_HPP_ @@ -26,25 +26,24 @@ #include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" #endif #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp" +#include "AcousticVTIFields.hpp" namespace geos { -using namespace fields; - /// Namespace to contain the acoustic wave kernels. -namespace acousticVTIFletcherAdjointWaveEquationSEMKernels +namespace acousticVTIFletcherWaveEquationSEMKernels { /** - * @brief Implements kernels for solving the acoustic wave equations + * @brief Implements kernels for solving the pseudo-acoustic VTI wave equations * explicit central FD method and SEM * @copydoc geos::finiteElement::KernelBase * @tparam SUBREGION_TYPE The type of subregion that the kernel will act on. * - * ### AcousticWaveEquationSEMKernel Description + * ### AcousticVTIFletcherWaveEquationSEMKernel Description * Implements the KernelBase interface functions required for solving - * the acoustic wave equations using the + * the VTI pseudo-acoustic wave Fletcher's set of equations using the * "finite element kernel application" functions such as * geos::finiteElement::RegionBasedKernelApplication. * @@ -56,11 +55,11 @@ namespace acousticVTIFletcherAdjointWaveEquationSEMKernels template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > -class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< SUBREGION_TYPE, - CONSTITUTIVE_TYPE, - FE_TYPE, - 1, - 1 > +class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 1, + 1 > { public: @@ -94,14 +93,14 @@ class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< * @param dt The time interval for the step. * elements to be processed during this kernel launch. */ - ExplicitAcousticVTIFletcherAdjointSEM( NodeManager & nodeManager, - EdgeManager const & edgeManager, - FaceManager const & faceManager, - localIndex const targetRegionIndex, - SUBREGION_TYPE const & elementSubRegion, - FE_TYPE const & finiteElementSpace, - CONSTITUTIVE_TYPE & inputConstitutiveType, - real64 const dt ): + ExplicitAcousticSEM( NodeManager & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + real64 const dt ): Base( elementSubRegion, finiteElementSpace, inputConstitutiveType ), @@ -125,7 +124,7 @@ class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< /** * @copydoc geos::finiteElement::KernelBase::StackVariables * - * ### ExplicitAcousticVTIFletcherAdjointSEM Description + * ### ExplicitAcousticSEM Description * Adds a stack arrays for the nodal force, primary displacement variable, etc. */ struct StackVariables : Base::StackVariables @@ -143,7 +142,7 @@ class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< real32 stiffnessVectorLocal_p[ numNodesPerElem ]{}; real32 stiffnessVectorLocal_q[ numNodesPerElem ]{}; real32 invDensity; - real32 vti_epsi; // (1 + 2*epsilon) + real32 vti_epsi; // (1 + 2*epsilon) real32 vti_2delta; // 2*delta real32 vti_f; // f parameter }; @@ -202,7 +201,7 @@ class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< /** * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel * - * ### ExplicitAcousticVTIFletcherAdjointSEM Description + * ### ExplicitAcousticSEM Description * Calculates stiffness vector * */ @@ -212,7 +211,7 @@ class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< localIndex const q, StackVariables & stack ) const { - // Pseudo Stiffness xy + // (A_xy nabla u)((A_xy nabla v) m_finiteElementSpace.template computeStiffnessxyTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_pp = val * stack.invDensity *(-stack.vti_epsi) * m_p_n[m_elemsToNodes( k, j )]; @@ -224,7 +223,7 @@ class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< stack.stiffnessVectorLocal_q[i] += localIncrement_qq; } ); - // Pseudo-Stiffness z + // (A_z nabla u)((A_z nabla v) m_finiteElementSpace.template computeStiffnesszTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_pp = val * stack.invDensity * (stack.vti_f-1) * m_p_n[m_elemsToNodes( k, j )]; @@ -270,12 +269,14 @@ class ExplicitAcousticVTIFletcherAdjointSEM : public finiteElement::KernelBase< real64 const m_dt; }; -/// The factory used to construct a ExplicitAcousticWaveEquation kernel. + + +/// The factory used to construct a ExplicitAcousticVTIFletcherAdjoint kernel. using ExplicitAcousticVTIFletcherAdjointSEMFactory = finiteElement::KernelFactory< ExplicitAcousticVTIFletcherAdjointSEM, - real64 >; + real64 >; -} // namespace acousticVTIFletcherAdjointWaveEquationSEMKernels +} // namespace acousticVTIFletcherWaveEquationSEMKernels } // namespace geos diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp index 273ef48c4d6..236b9a8be9d 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp @@ -29,10 +29,11 @@ #include "mesh/ElementType.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "events/EventManager.hpp" #include "physicsSolvers/wavePropagation/shared/WaveSolverUtils.hpp" #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp" #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp" +#include "events/EventManager.hpp" +// #include "AcousticPMLSEMKernel.hpp" // Not working with VTI #include "physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp" namespace geos @@ -42,7 +43,7 @@ using namespace dataRepository; using namespace fields; AcousticVTIFletcherWaveEquationSEM::AcousticVTIFletcherWaveEquationSEM( const std::string & name, - Group * const parent ): + Group * const parent ): WaveSolverBase( name, parent ) { @@ -92,10 +93,10 @@ void AcousticVTIFletcherWaveEquationSEM::registerDataOnMesh( Group & meshBodies acousticvtifields::AcousticLateralSurfaceNodeIndicator, acousticvtifields::AcousticBottomSurfaceNodeIndicator >( getName() ); - /// register PML auxiliary variables only when a PML is specified in the xml + /// PML is not supported if( m_usePML ) { - GEOS_ERROR( "This option is not supported yet" ); + GEOS_ERROR( "This option (PML) is not supported" ); } FaceManager & faceManager = mesh.getFaceManager(); @@ -111,9 +112,8 @@ void AcousticVTIFletcherWaveEquationSEM::registerDataOnMesh( Group & meshBodies subRegion.registerField< acousticfields::AcousticDensity >( getName() ); subRegion.registerField< acousticfields::PartialGradient >( getName() ); - subRegion.registerField< acousticvtifields::AcousticEpsilon >( getName() ); subRegion.registerField< acousticvtifields::AcousticDelta >( getName() ); - subRegion.registerField< acousticvtifields::AcousticSigma >( getName() ); + subRegion.registerField< acousticvtifields::AcousticEpsilon >( getName() ); } ); } ); @@ -128,7 +128,7 @@ void AcousticVTIFletcherWaveEquationSEM::postInputInitialization() } void AcousticVTIFletcherWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) + arrayView1d< string const > const & regionNames ) { GEOS_MARK_FUNCTION; @@ -153,18 +153,6 @@ void AcousticVTIFletcherWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLe receiverConstants.setValues< EXEC_POLICY >( -1 ); receiverIsLocal.zero(); - arrayView2d< real32 > const sourceValue = m_sourceValue.toView(); - real64 dt = 0; - 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() ) - { - dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() ); - } - } - mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, [&]( localIndex const, localIndex const er, localIndex const esr, @@ -172,7 +160,7 @@ void AcousticVTIFletcherWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLe CellElementSubRegion & elementSubRegion ) { GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, - getDataContext() << ": Invalid type of element, the acoustic solver is designed for hexahedral meshes only (C3D8), using the SEM formulation", + getDataContext() << ": Invalid type of element, the acoustic VTI Fletcher solver is designed for hexahedral meshes only (C3D8), using the SEM formulation", InputError ); arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); @@ -211,47 +199,42 @@ void AcousticVTIFletcherWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLe receiverCoordinates, receiverIsLocal, receiverNodeIds, - receiverConstants, - sourceValue, - dt, - m_timeSourceFrequency, - m_timeSourceDelay, - m_rickerOrder ); + receiverConstants ); } } ); -/* elementSubRegion.faceList().freeOnDevice(); + elementSubRegion.faceList().freeOnDevice(); baseMesh.getElemManager().getRegion( er ).getSubRegion< CellElementSubRegion >( esr ).nodeList().freeOnDevice(); elementSubRegion.getElementCenter().freeOnDevice(); elementSubRegion.ghostRank().freeOnDevice(); - elementSubRegion.localToGlobalMap().freeOnDevice();*/ + elementSubRegion.localToGlobalMap().freeOnDevice(); } ); -/* baseMesh.getNodeManager().localToGlobalMap().freeOnDevice(); - baseMesh.getNodeManager().elementList().toView().freeOnDevice(); - baseMesh.getFaceManager().nodeList().toView().freeOnDevice(); - baseMesh.getNodeManager().referencePosition().freeOnDevice(); - m_sourceCoordinates.freeOnDevice(); - m_receiverCoordinates.freeOnDevice(); - facesToNodes.freeOnDevice(); - nodesToElements.freeOnDevice();*/ + baseMesh.getNodeManager().localToGlobalMap().freeOnDevice(); + baseMesh.getNodeManager().elementList().toView().freeOnDevice(); + baseMesh.getFaceManager().nodeList().toView().freeOnDevice(); + baseMesh.getNodeManager().referencePosition().freeOnDevice(); + facesToNodes.freeOnDevice(); + nodesToElements.freeOnDevice(); } -void AcousticVTIFletcherWaveEquationSEM::addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ) +void AcousticVTIFletcherWaveEquationSEM::addSourceToRightHandSide( real64 const & time_n, arrayView1d< real32 > const rhs ) { arrayView2d< localIndex const > const sourceNodeIds = m_sourceNodeIds.toViewConst(); arrayView2d< real64 const > const sourceConstants = m_sourceConstants.toViewConst(); arrayView1d< localIndex const > const sourceIsAccessible = m_sourceIsAccessible.toViewConst(); - arrayView2d< real32 const > const sourceValue = m_sourceValue.toViewConst(); - - GEOS_THROW_IF( cycleNumber > sourceValue.size( 0 ), - getDataContext() << ": Too many steps compared to array size", - std::runtime_error ); + real32 const timeSourceFrequency = m_timeSourceFrequency; + real32 const timeSourceDelay = m_timeSourceDelay; + localIndex const rickerOrder = m_rickerOrder; + bool useSourceWaveletTables = m_useSourceWaveletTables; + arrayView1d< TableFunction::KernelWrapper const > const sourceWaveletTableWrappers = m_sourceWaveletTableWrappers.toViewConst(); forAll< EXEC_POLICY >( sourceConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const isrc ) { if( sourceIsAccessible[isrc] == 1 ) { + real64 const srcValue = + useSourceWaveletTables ? sourceWaveletTableWrappers[ isrc ].compute( &time_n ) : WaveSolverUtils::evaluateRicker( time_n, timeSourceFrequency, timeSourceDelay, rickerOrder ); for( localIndex inode = 0; inode < sourceConstants.size( 1 ); ++inode ) { - real32 const localIncrement = sourceConstants[isrc][inode] * sourceValue[cycleNumber][isrc]; + real32 const localIncrement = sourceConstants[isrc][inode] * srcValue; RAJA::atomicAdd< ATOMIC_POLICY >( &rhs[sourceNodeIds[isrc][inode]], localIncrement ); } } @@ -329,8 +312,9 @@ void AcousticVTIFletcherWaveEquationSEM::initializePostInitialConditionsPreSubGr arrayView1d< real32 const > const vti_delta = elementSubRegion.getField< acousticvtifields::AcousticDelta >(); arrayView1d< real32 const > const vti_sigma = elementSubRegion.getField< acousticvtifields::AcousticSigma >(); - /// Partial gradient if gradient has to be computed + /// Partial gradient if gradient as to be computed arrayView1d< real32 > grad = elementSubRegion.getField< acousticfields::PartialGradient >(); + grad.zero(); finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) @@ -347,32 +331,54 @@ void AcousticVTIFletcherWaveEquationSEM::initializePostInitialConditionsPreSubGr AcousticMatricesSEM::DampingMatrix< FE_TYPE > kernelD( finiteElement ); kernelD.template computeVTIFletcherDampingMatrices< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), - nodeCoords, - elemsToFaces, - facesToNodes, - facesDomainBoundaryIndicator, - freeSurfaceFaceIndicator, - lateralSurfaceFaceIndicator, - bottomSurfaceFaceIndicator, - velocity, - density, - vti_epsilon, - vti_delta, - vti_sigma, - damping_pp, - damping_pq, - damping_qp, - damping_qq ); + nodeCoords, + elemsToFaces, + facesToNodes, + facesDomainBoundaryIndicator, + freeSurfaceFaceIndicator, + lateralSurfaceFaceIndicator, + bottomSurfaceFaceIndicator, + velocity, + density, + vti_epsilon, + vti_delta, + vti_sigma, + damping_pp, + damping_pq, + damping_qp, + damping_qq ); + } ); } ); } ); + if( m_timestepStabilityLimit==1 ) + { + // TODO: adapt to VTI + GEOS_ERROR( "This option (Time Step computation) is not supported" ); +/* real64 dtOut = 0.0; + computeTimeStep( dtOut ); + m_timestepStabilityLimit = 0; + m_timeStep=dtOut;*/ + + } + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); m_seismoCoeff.resize( m_receiverIsLocal.size()); m_seismoCoeff.setValues< EXEC_POLICY >( 0.5 ); } +//This function is only to give an easy accesss to the computation of the time-step for Pygeosx interface and avoid to exit the code when +// using Pygeosx + +real64 AcousticVTIFletcherWaveEquationSEM::computeTimeStep( real64 & dtOut ) +{ + // TODO: adapt to VTI + GEOS_ERROR( "This option (Time Step computation) is not supported" ); +} + + void AcousticVTIFletcherWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain ) { real64 const time = 0.0; @@ -385,14 +391,14 @@ void AcousticVTIFletcherWaveEquationSEM::precomputeSurfaceFieldIndicator( Domain ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); /// array of indicators: 1 if a face is on on lateral surface; 0 otherwise - arrayView1d< localIndex > const lateralSurfaceFaceIndicator = faceManager.getField< fields::acousticvtifields::AcousticLateralSurfaceFaceIndicator >(); + arrayView1d< localIndex > const lateralSurfaceFaceIndicator = faceManager.getField< acoustivtifields::AcousticLateralSurfaceFaceIndicator >(); /// array of indicators: 1 if a node is on on lateral surface; 0 otherwise - arrayView1d< localIndex > const lateralSurfaceNodeIndicator = nodeManager.getField< fields::acousticvtifields::AcousticLateralSurfaceNodeIndicator >(); + arrayView1d< localIndex > const lateralSurfaceNodeIndicator = nodeManager.getField< acoustivtifields::AcousticLateralSurfaceNodeIndicator >(); /// array of indicators: 1 if a face is on on bottom surface; 0 otherwise - arrayView1d< localIndex > const bottomSurfaceFaceIndicator = faceManager.getField< fields::acousticvtifields::AcousticBottomSurfaceFaceIndicator >(); + arrayView1d< localIndex > const bottomSurfaceFaceIndicator = faceManager.getField< acoustivtifields::AcousticBottomSurfaceFaceIndicator >(); /// array of indicators: 1 if a node is on on bottom surface; 0 otherwise - arrayView1d< localIndex > const bottomSurfaceNodeIndicator = nodeManager.getField< fields::acousticvtifields::AcousticBottomSurfaceNodeIndicator >(); + arrayView1d< localIndex > const bottomSurfaceNodeIndicator = nodeManager.getField< acoustivtifields::AcousticBottomSurfaceNodeIndicator >(); // Lateral surfaces fsManager.apply< FaceManager >( time, @@ -486,9 +492,6 @@ void AcousticVTIFletcherWaveEquationSEM::applyFreeSurfaceBC( real64 time, Domain /// array of indicators: 1 if a node is on on free surface; 0 otherwise arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfields::AcousticFreeSurfaceNodeIndicator >(); - // freeSurfaceFaceIndicator.zero(); - // freeSurfaceNodeIndicator.zero(); - fsManager.apply< FaceManager >( time, domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), string( "FreeSurface" ), @@ -533,26 +536,119 @@ void AcousticVTIFletcherWaveEquationSEM::applyFreeSurfaceBC( real64 time, Domain } void AcousticVTIFletcherWaveEquationSEM::initializePML() -{ - GEOS_ERROR( "This option is not supported yet" ); +{ GEOS_ERROR( "This option (PML) is not supported" ); return; } -void AcousticVTIFletcherWaveEquationSEM::applyPML( real64 const GEOS_UNUSED_PARAM( time ), DomainPartition & GEOS_UNUSED_PARAM( domain ) ) +void AcousticVTIFletcherWaveEquationSEM::applyPML( real64 const time, DomainPartition & domain ) { - GEOS_ERROR( "This option is not supported yet" ); - return; + GEOS_MARK_FUNCTION; + + FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); + parametersPML const & param = getReference< parametersPML >( viewKeyStruct::parametersPMLString() ); + + /// Loop over the different mesh bodies; for wave propagation, there is only one mesh body + /// which is the whole mesh + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & ) + { + + NodeManager & nodeManager = mesh.getNodeManager(); + + /// Array views of the pressure p, PML auxiliary variables, and node coordinates + arrayView1d< real32 const > const p_n = nodeManager.getField< acousticfields::Pressure_n >(); + arrayView2d< real32 const > const v_n = nodeManager.getField< acousticfields::AuxiliaryVar1PML >(); + arrayView2d< real32 > const grad_n = nodeManager.getField< acousticfields::AuxiliaryVar2PML >(); + arrayView1d< real32 > const divV_n = nodeManager.getField< acousticfields::AuxiliaryVar3PML >(); + arrayView1d< real32 const > const u_n = nodeManager.getField< acousticfields::AuxiliaryVar4PML >(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords32 = nodeManager.getField< fields::referencePosition32 >().toViewConst(); + + /// Select the subregions concerned by the PML (specified in the xml by the Field Specification) + /// 'targetSet' contains the indices of the elements in a given subregion + fsManager.apply< ElementSubRegionBase, + PerfectlyMatchedLayer >( time, + mesh, + PerfectlyMatchedLayer::catalogName(), + [&]( PerfectlyMatchedLayer const &, + string const &, + SortedArrayView< localIndex const > const & targetSet, + ElementSubRegionBase & subRegion, + string const & ) + + { + + /// Get the element to nodes mapping in the subregion + CellElementSubRegion::NodeMapType const & elemToNodes = + subRegion.getReference< CellElementSubRegion::NodeMapType >( CellElementSubRegion::viewKeyStruct::nodeListString() ); + + /// Get a const ArrayView of the mapping above + traits::ViewTypeConst< CellElementSubRegion::NodeMapType > const elemToNodesViewConst = elemToNodes.toViewConst(); + + /// Array view of the wave speed + arrayView1d< real32 const > const vel = subRegion.getReference< array1d< real32 > >( acousticfields::AcousticVelocity::key()); + + /// Get the object needed to determine the type of the element in the subregion + finiteElement::FiniteElementBase const & + fe = subRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + + real32 xMin[3]; + real32 xMax[3]; + real32 dMin[3]; + real32 dMax[3]; + real32 cMin[3]; + real32 cMax[3]; + for( integer i=0; i<3; ++i ) + { + xMin[i] = param.xMinPML[i]; + xMax[i] = param.xMaxPML[i]; + dMin[i] = param.thicknessMinXYZPML[i]; + dMax[i] = param.thicknessMaxXYZPML[i]; + cMin[i] = param.waveSpeedMinXYZPML[i]; + cMax[i] = param.waveSpeedMaxXYZPML[i]; + } + real32 const r = param.reflectivityPML; + + /// Get the type of the elements in the subregion + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) + { + using FE_TYPE = TYPEOFREF( finiteElement ); + + /// apply the PML kernel + AcousticPMLSEM:: + PMLKernel< FE_TYPE > kernel( finiteElement ); + kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > + ( targetSet, + nodeCoords32, + elemToNodesViewConst, + vel, + p_n, + v_n, + u_n, + xMin, + xMax, + dMin, + dMax, + cMin, + cMax, + r, + grad_n, + divV_n ); + } ); + } ); + } ); + } real64 AcousticVTIFletcherWaveEquationSEM::explicitStepForward( real64 const & time_n, - real64 const & dt, - integer cycleNumber, - DomainPartition & domain, - bool computeGradient ) + real64 const & dt, + integer cycleNumber, + DomainPartition & domain, + bool computeGradient ) { - real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain, true ); + real64 dtCompute = explicitStepInternal( time_n, dt, domain, true ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -577,26 +673,26 @@ real64 AcousticVTIFletcherWaveEquationSEM::explicitStepForward( real64 const & t prepareNextTimestep( mesh ); } ); - return dtOut; + return dtCompute; } real64 AcousticVTIFletcherWaveEquationSEM::explicitStepBackward( real64 const & time_n, - real64 const & dt, - integer cycleNumber, - DomainPartition & domain, - bool computeGradient ) + real64 const & dt, + integer cycleNumber, + DomainPartition & domain, + bool computeGradient ) { - - real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain, false ); - + real64 dtCompute = explicitStepInternal( time_n, dt, domain, false ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, - arrayView1d< string const > const & GEOS_UNUSED_PARAM ( regionNames ) ) + arrayView1d< string const > const & regionNames ) { NodeManager & nodeManager = mesh.getNodeManager(); + arrayView1d< real32 const > const mass = nodeManager.getField< acousticfields::AcousticMassVector >(); + arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticvtifields::Pressure_p_nm1 >(); arrayView1d< real32 > const p_n = nodeManager.getField< acousticvtifields::Pressure_p_n >(); arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticvtifields::Pressure_p_np1 >(); @@ -605,6 +701,11 @@ real64 AcousticVTIFletcherWaveEquationSEM::explicitStepBackward( real64 const & arrayView1d< real32 > const q_n = nodeManager.getField< acousticvtifields::Pressure_q_n >(); arrayView1d< real32 > const q_np1 = nodeManager.getField< acousticvtifields::Pressure_q_np1 >(); + + EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); + real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() ); + int const maxCycle = int(round( maxTime / dt )); + if( computeGradient && cycleNumber >= 0 ) { GEOS_ERROR( "This option is not supported yet" ); @@ -613,7 +714,7 @@ real64 AcousticVTIFletcherWaveEquationSEM::explicitStepBackward( real64 const & prepareNextTimestep( mesh ); } ); - return dtOut; + return dtCompute; } void AcousticVTIFletcherWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh ) @@ -649,13 +750,12 @@ void AcousticVTIFletcherWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh ) } ); } -void AcousticVTIFletcherWaveEquationSEM::computeUnknowns( real64 const & GEOS_UNUSED_PARAM( time_n ), - real64 const & dt, - integer cycleNumber, - DomainPartition & GEOS_UNUSED_PARAM( domain ), - MeshLevel & mesh, - arrayView1d< string const > const & regionNames, - bool const isForward ) +void AcousticVTIFletcherWaveEquationSEM::computeUnknowns( real64 const & time_n, + real64 const & dt, + DomainPartition & domain, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames, + bool const isForward ) { NodeManager & nodeManager = mesh.getNodeManager(); @@ -673,13 +773,17 @@ void AcousticVTIFletcherWaveEquationSEM::computeUnknowns( real64 const & GEOS_UN arrayView1d< real32 > const q_n = nodeManager.getField< acousticvtifields::Pressure_q_n >(); arrayView1d< real32 > const q_np1 = nodeManager.getField< acousticvtifields::Pressure_q_np1 >(); + arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticfields::Pressure_nm1 >(); + arrayView1d< real32 > const p_n = nodeManager.getField< acousticfields::Pressure_n >(); + arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticfields::Pressure_np1 >(); + arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfields::AcousticFreeSurfaceNodeIndicator >(); arrayView1d< localIndex const > const lateralSurfaceNodeIndicator = nodeManager.getField< acousticvtifields::AcousticLateralSurfaceNodeIndicator >(); arrayView1d< localIndex const > const bottomSurfaceNodeIndicator = nodeManager.getField< acousticvtifields::AcousticBottomSurfaceNodeIndicator >(); arrayView1d< real32 > const stiffnessVector_p = nodeManager.getField< acousticvtifields::StiffnessVector_p >(); arrayView1d< real32 > const stiffnessVector_q = nodeManager.getField< acousticvtifields::StiffnessVector_q >(); arrayView1d< real32 > const rhs = nodeManager.getField< acousticfields::ForcingRHS >(); - if( isForward ) +if( isForward ) { auto kernelFactory = acousticVTIFletcherWaveEquationSEMKernels::ExplicitAcousticVTIFletcherSEMFactory( dt ); @@ -704,22 +808,20 @@ void AcousticVTIFletcherWaveEquationSEM::computeUnknowns( real64 const & GEOS_UN getDiscretizationName(), "", kernelFactory ); - } + } //Modification of cycleNember useful when minTime < 0 - EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); - real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() ); - integer const cycleForSource = int(round( -minTime / dt + cycleNumber )); - addSourceToRightHandSide( cycleForSource, rhs ); + addSourceToRightHandSide( time_n, rhs ); + + /// calculate your time integrators + real64 const dt2 = pow( dt, 2 ); SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst(); if( !m_usePML ) { GEOS_MARK_SCOPE ( updateP ); - AcousticTimeSchemeSEM::LeapFrogVTIWithoutPML( dt, p_np1, p_n, p_nm1, q_np1, q_n, q_nm1, mass, stiffnessVector_p, - stiffnessVector_q, damping_pp, damping_pq, damping_qp, damping_qq, - rhs, freeSurfaceNodeIndicator, lateralSurfaceNodeIndicator, - bottomSurfaceNodeIndicator, solverTargetNodesSet ); + AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, + rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); } else { @@ -728,11 +830,10 @@ void AcousticVTIFletcherWaveEquationSEM::computeUnknowns( real64 const & GEOS_UN } void AcousticVTIFletcherWaveEquationSEM::synchronizeUnknowns( real64 const & time_n, - real64 const & dt, - integer const, - DomainPartition & domain, - MeshLevel & mesh, - arrayView1d< string const > const & ) + real64 const & dt, + DomainPartition & domain, + MeshLevel & mesh, + arrayView1d< string const > const & ) { NodeManager & nodeManager = mesh.getNodeManager(); @@ -755,6 +856,7 @@ void AcousticVTIFletcherWaveEquationSEM::synchronizeUnknowns( real64 const & tim GEOS_ERROR( "This option is not supported yet" ); } + CommunicationTools & syncFields = CommunicationTools::getInstance(); syncFields.synchronizeFields( fieldsToBeSync, mesh, @@ -763,7 +865,6 @@ void AcousticVTIFletcherWaveEquationSEM::synchronizeUnknowns( real64 const & tim /// compute the seismic traces since last step. arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); - // Seismo trace = (p+q)/2 computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers, m_seismoCoeff.toView(), false ); computeAllSeismoTraces( time_n, dt, q_np1, q_n, pReceivers, m_seismoCoeff.toView(), true ); incrementIndexSeismoTrace( time_n ); @@ -775,34 +876,36 @@ void AcousticVTIFletcherWaveEquationSEM::synchronizeUnknowns( real64 const & tim } real64 AcousticVTIFletcherWaveEquationSEM::explicitStepInternal( real64 const & time_n, - real64 const & dt, - integer const cycleNumber, - DomainPartition & domain, - bool const isForward ) + real64 const & dt, + DomainPartition & domain, + bool const isForward ) { GEOS_MARK_FUNCTION; GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc ); + real64 dtCompute; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - computeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames, isForward ); - synchronizeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames ); + localIndex nSubSteps = (int) ceil( dt/m_timeStep ); + dtCompute = dt/nSubSteps; + computeUnknowns( time_n, dtCompute, domain, mesh, regionNames, isForward ); + synchronizeUnknowns( time_n, dtCompute, domain, mesh, regionNames ); } ); - return dt; + return dtCompute; } void AcousticVTIFletcherWaveEquationSEM::cleanup( real64 const time_n, - integer const cycleNumber, - integer const eventCounter, - real64 const eventProgress, - DomainPartition & domain ) + integer const cycleNumber, + integer const eventCounter, + real64 const eventProgress, + DomainPartition & domain ) { // call the base class cleanup (for reporting purposes) - SolverBase::cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain ); + PhysicsSolverBase::cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain ); // compute the remaining seismic traces, if needed forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -823,6 +926,6 @@ void AcousticVTIFletcherWaveEquationSEM::cleanup( real64 const time_n, } ); } -REGISTER_CATALOG_ENTRY( SolverBase, AcousticVTIFletcherWaveEquationSEM, string const &, dataRepository::Group * const ) +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, AcousticVTIFletcherWaveEquationSEM, string const &, dataRepository::Group * const ) } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.hpp index 69066110960..37a531b90cc 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.hpp @@ -23,7 +23,7 @@ #include "physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp" #include "mesh/MeshFields.hpp" -#include "physicsSolvers/SolverBase.hpp" +#include "physicsSolvers/PhysicsSolverBase.hpp" #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp" #include "AcousticVTIFields.hpp" @@ -38,7 +38,7 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase using ATOMIC_POLICY = AtomicPolicy< EXEC_POLICY >; AcousticVTIFletcherWaveEquationSEM( const std::string & name, - Group * const parent ); + Group * const parent ); virtual ~AcousticVTIFletcherWaveEquationSEM() override; @@ -54,7 +54,7 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase static string catalogName() { return "AcousticVTIFletcherSEM"; } /** - * @copydoc SolverBase::getCatalogName() + * @copydoc PhysicsSolverBase::getCatalogName() */ string getCatalogName() const override { return catalogName(); } @@ -88,7 +88,7 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase * @param cycleNumber the cycle number/step number of evaluation of the source * @param rhs the right hand side vector to be computed */ - virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ); + virtual void addSourceToRightHandSide( real64 const & time_n, arrayView1d< real32 > const rhs ); /** @@ -96,6 +96,11 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase */ virtual void initializePML() override; + /** + */ + virtual real64 computeTimeStep( real64 & dtOut ) override; + + /** * @brief Overridden from ExecutableGroup. Used to write last seismogram if needed. @@ -119,13 +124,11 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase */ real64 explicitStepInternal( real64 const & time_n, real64 const & dt, - integer const cycleNumber, DomainPartition & domain, bool const isForward ); void computeUnknowns( real64 const & time_n, real64 const & dt, - integer const cycleNumber, DomainPartition & domain, MeshLevel & mesh, arrayView1d< string const > const & regionNames, @@ -133,7 +136,6 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase void synchronizeUnknowns( real64 const & time_n, real64 const & dt, - integer const cycleNumber, DomainPartition & domain, MeshLevel & mesh, arrayView1d< string const > const & regionNames ); @@ -156,12 +158,6 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase */ virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - /** - * @brief Compute the lateral and bottom surface Field indicators of the boxed domain - * @param domain the partition domain - */ - virtual void precomputeSurfaceFieldIndicator( DomainPartition & domain ); - /** * @brief Apply free surface condition to the face define in the geometry box from the xml * @param time the time to apply the BC @@ -179,7 +175,7 @@ class AcousticVTIFletcherWaveEquationSEM : public WaveSolverBase /// Pressure_np1 at the receiver location for each time step for each receiver array2d< real32 > m_pressureNp1AtReceivers; - /// Array of size the number of receivers and full of 0.5 (used for calculating the seismos) + /// Array of size the number of receivers and filled with 0.5 (used for calculating the seismos) array1d< real32 > m_seismoCoeff; }; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEMKernel.hpp index c5ea452a771..3ae04ae309a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEMKernel.hpp @@ -26,18 +26,17 @@ #include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" #endif #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp" +#include "AcousticVTIFields.hpp" namespace geos { -using namespace fields; - /// Namespace to contain the acoustic wave kernels. namespace acousticVTIFletcherWaveEquationSEMKernels { /** - * @brief Implements kernels for solving the acoustic wave equations + * @brief Implements kernels for solving the pseudo-acoustic VTI wave equations * explicit central FD method and SEM * @copydoc geos::finiteElement::KernelBase * @tparam SUBREGION_TYPE The type of subregion that the kernel will act on. @@ -56,11 +55,11 @@ namespace acousticVTIFletcherWaveEquationSEMKernels template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > -class ExplicitAcousticVTIFletcherSEM : public finiteElement::KernelBase< SUBREGION_TYPE, - CONSTITUTIVE_TYPE, - FE_TYPE, - 1, - 1 > +class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 1, + 1 > { public: @@ -94,14 +93,14 @@ class ExplicitAcousticVTIFletcherSEM : public finiteElement::KernelBase< SUBREGI * @param dt The time interval for the step. * elements to be processed during this kernel launch. */ - ExplicitAcousticVTIFletcherSEM( NodeManager & nodeManager, - EdgeManager const & edgeManager, - FaceManager const & faceManager, - localIndex const targetRegionIndex, - SUBREGION_TYPE const & elementSubRegion, - FE_TYPE const & finiteElementSpace, - CONSTITUTIVE_TYPE & inputConstitutiveType, - real64 const dt ): + ExplicitAcousticSEM( NodeManager & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + real64 const dt ): Base( elementSubRegion, finiteElementSpace, inputConstitutiveType ), @@ -125,7 +124,7 @@ class ExplicitAcousticVTIFletcherSEM : public finiteElement::KernelBase< SUBREGI /** * @copydoc geos::finiteElement::KernelBase::StackVariables * - * ### ExplicitAcousticVTIFletcherSEM Description + * ### ExplicitAcousticSEM Description * Adds a stack arrays for the nodal force, primary displacement variable, etc. */ struct StackVariables : Base::StackVariables @@ -202,7 +201,7 @@ class ExplicitAcousticVTIFletcherSEM : public finiteElement::KernelBase< SUBREGI /** * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel * - * ### ExplicitAcousticVTIFletcherSEM Description + * ### ExplicitAcousticSEM Description * Calculates stiffness vector * */ @@ -212,17 +211,17 @@ class ExplicitAcousticVTIFletcherSEM : public finiteElement::KernelBase< SUBREGI localIndex const q, StackVariables & stack ) const { - // Pseudo Stiffness xy + // (A_xy nabla u)((A_xy nabla v) m_finiteElementSpace.template computeStiffnessxyTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_p = val * stack.invDensity *(-stack.vti_epsi) * m_p_n[m_elemsToNodes( k, j )]; stack.stiffnessVectorLocal_p[i] += localIncrement_p; + real32 const localIncrement_q = val * stack.invDensity * ((-stack.vti_2delta-stack.vti_f) * m_p_n[m_elemsToNodes( k, j )] +(stack.vti_f-1)*m_q_n[m_elemsToNodes( k, j )]); stack.stiffnessVectorLocal_q[i] += localIncrement_q; } ); - // Pseudo-Stiffness z - + // (A_z nabla u)((A_z nabla v) m_finiteElementSpace.template computeStiffnesszTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_p = val * stack.invDensity * ((stack.vti_f-1)*m_p_n[m_elemsToNodes( k, j )] - stack.vti_f*m_q_n[m_elemsToNodes( k, j )]); @@ -267,9 +266,9 @@ class ExplicitAcousticVTIFletcherSEM : public finiteElement::KernelBase< SUBREGI -/// The factory used to construct a ExplicitAcousticWaveEquation kernel. +/// The factory used to construct a ExplicitAcousticVTIFletcher kernel. using ExplicitAcousticVTIFletcherSEMFactory = finiteElement::KernelFactory< ExplicitAcousticVTIFletcherSEM, - real64 >; + real64 >; } // namespace acousticVTIFletcherWaveEquationSEMKernels diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangAdjointWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangAdjointWaveEquationSEMKernel.hpp index d57563570e3..08c5c0d9f15 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangAdjointWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangAdjointWaveEquationSEMKernel.hpp @@ -26,25 +26,24 @@ #include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" #endif #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp" +#include "AcousticVTIFields.hpp" namespace geos { -using namespace fields; - /// Namespace to contain the acoustic wave kernels. -namespace acousticVTIZhangAdjointWaveEquationSEMKernels +namespace acousticVTIZhangWaveEquationSEMKernels { /** - * @brief Implements kernels for solving the acoustic wave equations + * @brief Implements kernels for solving the pseudo-acoustic VTI wave equations * explicit central FD method and SEM * @copydoc geos::finiteElement::KernelBase * @tparam SUBREGION_TYPE The type of subregion that the kernel will act on. * - * ### AcousticVTIZhangAdjointWaveEquationSEMKernels Description + * ### AcousticVTIZhangWaveEquationSEMKernel Description * Implements the KernelBase interface functions required for solving - * the Adjoint of the VTI pseudo-acoustic wave Zhang's set of equations using the + * the VTI pseudo-acoustic wave Zhang's set of equations using the * "finite element kernel application" functions such as * geos::finiteElement::RegionBasedKernelApplication. * @@ -56,11 +55,11 @@ namespace acousticVTIZhangAdjointWaveEquationSEMKernels template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > -class ExplicitAcousticVTIZhangAdjointSEM : public finiteElement::KernelBase< SUBREGION_TYPE, - CONSTITUTIVE_TYPE, - FE_TYPE, - 1, - 1 > +class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 1, + 1 > { public: @@ -94,14 +93,14 @@ class ExplicitAcousticVTIZhangAdjointSEM : public finiteElement::KernelBase< SUB * @param dt The time interval for the step. * elements to be processed during this kernel launch. */ - ExplicitAcousticVTIZhangAdjointSEM( NodeManager & nodeManager, - EdgeManager const & edgeManager, - FaceManager const & faceManager, - localIndex const targetRegionIndex, - SUBREGION_TYPE const & elementSubRegion, - FE_TYPE const & finiteElementSpace, - CONSTITUTIVE_TYPE & inputConstitutiveType, - real64 const dt ): + ExplicitAcousticSEM( NodeManager & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + real64 const dt ): Base( elementSubRegion, finiteElementSpace, inputConstitutiveType ), @@ -124,7 +123,7 @@ class ExplicitAcousticVTIZhangAdjointSEM : public finiteElement::KernelBase< SUB /** * @copydoc geos::finiteElement::KernelBase::StackVariables * - * ### ExplicitAcousticVTIZhangAdjointSEM Description + * ### ExplicitAcousticSEM Description * Adds a stack arrays for the nodal force, primary displacement variable, etc. */ struct StackVariables : Base::StackVariables @@ -199,7 +198,7 @@ class ExplicitAcousticVTIZhangAdjointSEM : public finiteElement::KernelBase< SUB /** * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel * - * ### ExplicitAcousticVTIZhangAdjointSEM Description + * ### ExplicitAcousticSEM Description * Calculates stiffness vector * */ @@ -209,7 +208,7 @@ class ExplicitAcousticVTIZhangAdjointSEM : public finiteElement::KernelBase< SUB localIndex const q, StackVariables & stack ) const { - // Pseudo Stiffness xy + // (A_xy nabla u)((A_xy nabla v) m_finiteElementSpace.template computeStiffnessxyTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_pp = -val * stack.invDensity * stack.vti_epsi * m_p_n[m_elemsToNodes( k, j )]; @@ -219,7 +218,7 @@ class ExplicitAcousticVTIZhangAdjointSEM : public finiteElement::KernelBase< SUB stack.stiffnessVectorLocal_p[i] += localIncrement_pq; } ); - // Pseudo-Stiffness z + // (A_z nabla u)((A_z nabla v) m_finiteElementSpace.template computeStiffnesszTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_qp = -val * stack.invDensity * stack.vti_sqrtDelta * m_p_n[m_elemsToNodes( k, j )]; @@ -261,12 +260,12 @@ class ExplicitAcousticVTIZhangAdjointSEM : public finiteElement::KernelBase< SUB -/// The factory used to construct a ExplicitAcousticWaveEquation kernel. +/// The factory used to construct a ExplicitAcousticVTIZhangAdjoint kernel. using ExplicitAcousticVTIZhangAdjointSEMFactory = finiteElement::KernelFactory< ExplicitAcousticVTIZhangAdjointSEM, real64 >; -} // namespace acousticVTIZhangAdjointWaveEquationSEMKernels +} // namespace acousticVTIZhangWaveEquationSEMKernels } // namespace geos diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp index 442653e5a4e..3879e3e6748 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp @@ -29,10 +29,11 @@ #include "mesh/ElementType.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "events/EventManager.hpp" #include "physicsSolvers/wavePropagation/shared/WaveSolverUtils.hpp" #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticTimeSchemeSEMKernel.hpp" #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp" +#include "events/EventManager.hpp" +// #include "AcousticPMLSEMKernel.hpp" // Not working with VTI #include "physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp" namespace geos @@ -42,7 +43,7 @@ using namespace dataRepository; using namespace fields; AcousticVTIZhangWaveEquationSEM::AcousticVTIZhangWaveEquationSEM( const std::string & name, - Group * const parent ): + Group * const parent ): WaveSolverBase( name, parent ) { @@ -92,10 +93,10 @@ void AcousticVTIZhangWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) acousticvtifields::AcousticLateralSurfaceNodeIndicator, acousticvtifields::AcousticBottomSurfaceNodeIndicator >( getName() ); - /// register PML auxiliary variables only when a PML is specified in the xml + /// PML is not supported if( m_usePML ) { - GEOS_ERROR( "This option is not supported yet" ); + GEOS_ERROR( "This option (PML) is not supported" ); } FaceManager & faceManager = mesh.getFaceManager(); @@ -127,7 +128,7 @@ void AcousticVTIZhangWaveEquationSEM::postInputInitialization() } void AcousticVTIZhangWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) + arrayView1d< string const > const & regionNames ) { GEOS_MARK_FUNCTION; @@ -159,7 +160,7 @@ void AcousticVTIZhangWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel CellElementSubRegion & elementSubRegion ) { GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, - getDataContext() << ": Invalid type of element, the acoustic solver is designed for hexahedral meshes only (C3D8), using the SEM formulation", + getDataContext() << ": Invalid type of element, the acoustic VTI Zhang solver is designed for hexahedral meshes only (C3D8), using the SEM formulation", InputError ); arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); @@ -175,7 +176,6 @@ void AcousticVTIZhangWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel { using FE_TYPE = TYPEOFREF( finiteElement ); -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp { GEOS_MARK_SCOPE( acousticVTIZhangWaveEquationSEMKernels::PrecomputeSourceAndReceiverKernel ); PreComputeSourcesAndReceivers:: @@ -199,57 +199,15 @@ void AcousticVTIZhangWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel receiverCoordinates, receiverIsLocal, receiverNodeIds, - receiverConstants, - sourceValue, - dt, - m_timeSourceFrequency, - m_timeSourceDelay, - m_rickerOrder ); + receiverConstants ); } -======= - acousticVTIWaveEquationSEMKernels:: - PrecomputeSourceAndReceiverKernel:: - launch< EXEC_POLICY, FE_TYPE > - ( elementSubRegion.size(), - facesToNodes, - nodeCoords, - nodeLocalToGlobal, - elemLocalToGlobal, - nodesToElements, - baseElemsToNodes, - elemGhostRank, - elemsToNodes, - elemsToFaces, - elemCenter, - sourceCoordinates, - sourceIsAccessible, - sourceNodeIds, - sourceConstants, - receiverCoordinates, - receiverIsLocal, - receiverNodeIds, - receiverConstants ); ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp } ); -/* elementSubRegion.faceList().freeOnDevice(); + elementSubRegion.faceList().freeOnDevice(); baseMesh.getElemManager().getRegion( er ).getSubRegion< CellElementSubRegion >( esr ).nodeList().freeOnDevice(); elementSubRegion.getElementCenter().freeOnDevice(); elementSubRegion.ghostRank().freeOnDevice(); - elementSubRegion.localToGlobalMap().freeOnDevice();*/ + elementSubRegion.localToGlobalMap().freeOnDevice(); } ); -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp -/* baseMesh.getNodeManager().localToGlobalMap().freeOnDevice(); - baseMesh.getNodeManager().elementList().toView().freeOnDevice(); - baseMesh.getFaceManager().nodeList().toView().freeOnDevice(); - baseMesh.getNodeManager().referencePosition().freeOnDevice(); - m_sourceCoordinates.freeOnDevice(); - m_receiverCoordinates.freeOnDevice(); - facesToNodes.freeOnDevice(); - nodesToElements.freeOnDevice();*/ -} - -void AcousticVTIZhangWaveEquationSEM::addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ) -======= baseMesh.getNodeManager().localToGlobalMap().freeOnDevice(); baseMesh.getNodeManager().elementList().toView().freeOnDevice(); baseMesh.getFaceManager().nodeList().toView().freeOnDevice(); @@ -258,28 +216,22 @@ void AcousticVTIZhangWaveEquationSEM::addSourceToRightHandSide( integer const & nodesToElements.freeOnDevice(); } -void AcousticVTIWaveEquationSEM::addSourceToRightHandSide( real64 const & time_n, arrayView1d< real32 > const rhs ) ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +void AcousticVTIZhangWaveEquationSEM::addSourceToRightHandSide( real64 const & time_n, arrayView1d< real32 > const rhs ) { arrayView2d< localIndex const > const sourceNodeIds = m_sourceNodeIds.toViewConst(); arrayView2d< real64 const > const sourceConstants = m_sourceConstants.toViewConst(); arrayView1d< localIndex const > const sourceIsAccessible = m_sourceIsAccessible.toViewConst(); -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp - arrayView2d< real32 const > const sourceValue = m_sourceValue.toViewConst(); - - GEOS_THROW_IF( cycleNumber > sourceValue.size( 0 ), - getDataContext() << ": Too many steps compared to array size", - std::runtime_error ); -======= + real32 const timeSourceFrequency = m_timeSourceFrequency; + real32 const timeSourceDelay = m_timeSourceDelay; + localIndex const rickerOrder = m_rickerOrder; bool useSourceWaveletTables = m_useSourceWaveletTables; - real64 const rickerValue = useSourceWaveletTables ? 0 : WaveSolverUtils::evaluateRicker( time_n, m_timeSourceFrequency, m_timeSourceDelay, m_rickerOrder ); arrayView1d< TableFunction::KernelWrapper const > const sourceWaveletTableWrappers = m_sourceWaveletTableWrappers.toViewConst(); ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp forAll< EXEC_POLICY >( sourceConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const isrc ) { if( sourceIsAccessible[isrc] == 1 ) { - real64 const srcValue = useSourceWaveletTables ? sourceWaveletTableWrappers[ isrc ].compute( &time_n ) : rickerValue; + real64 const srcValue = + useSourceWaveletTables ? sourceWaveletTableWrappers[ isrc ].compute( &time_n ) : WaveSolverUtils::evaluateRicker( time_n, timeSourceFrequency, timeSourceDelay, rickerOrder ); for( localIndex inode = 0; inode < sourceConstants.size( 1 ); ++inode ) { real32 const localIncrement = sourceConstants[isrc][inode] * srcValue; @@ -359,8 +311,9 @@ void AcousticVTIZhangWaveEquationSEM::initializePostInitialConditionsPreSubGroup arrayView1d< real32 const > const vti_epsilon = elementSubRegion.getField< acousticvtifields::AcousticEpsilon >(); arrayView1d< real32 const > const vti_delta = elementSubRegion.getField< acousticvtifields::AcousticDelta >(); - /// Partial gradient if gradient has to be computed + /// Partial gradient if gradient as to be computed arrayView1d< real32 > grad = elementSubRegion.getField< acousticfields::PartialGradient >(); + grad.zero(); finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) @@ -393,27 +346,38 @@ void AcousticVTIZhangWaveEquationSEM::initializePostInitialConditionsPreSubGroup damping_qp, damping_qq ); + } ); } ); } ); + if( m_timestepStabilityLimit==1 ) + { + // TODO: adapt to VTI + GEOS_ERROR( "This option (Time Step computation) is not supported" ); +/* real64 dtOut = 0.0; + computeTimeStep( dtOut ); + m_timestepStabilityLimit = 0; + m_timeStep=dtOut;*/ + + } + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal ); m_seismoCoeff.resize( m_receiverIsLocal.size()); m_seismoCoeff.setValues< EXEC_POLICY >( 0.5 ); } -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp -void AcousticVTIZhangWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain ) -======= -real64 AcousticVTIWaveEquationSEM::computeTimeStep( real64 & dtOut ) +//This function is only to give an easy accesss to the computation of the time-step for Pygeosx interface and avoid to exit the code when +// using Pygeosx + +real64 AcousticVTIZhangWaveEquationSEM::computeTimeStep( real64 & dtOut ) { - GEOS_ERROR( getDataContext() << ": Time-Step computation for the second order acoustic vti wave propagator not yet implemented" ); - return dtOut; + // TODO: adapt to VTI + GEOS_ERROR( "This option (Time Step computation) is not supported" ); } -void AcousticVTIWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain ) ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +void AcousticVTIZhangWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain ) { real64 const time = 0.0; FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); @@ -425,14 +389,14 @@ void AcousticVTIWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartitio ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); /// array of indicators: 1 if a face is on on lateral surface; 0 otherwise - arrayView1d< localIndex > const lateralSurfaceFaceIndicator = faceManager.getField< fields::acousticvtifields::AcousticLateralSurfaceFaceIndicator >(); + arrayView1d< localIndex > const lateralSurfaceFaceIndicator = faceManager.getField< acoustivtifields::AcousticLateralSurfaceFaceIndicator >(); /// array of indicators: 1 if a node is on on lateral surface; 0 otherwise - arrayView1d< localIndex > const lateralSurfaceNodeIndicator = nodeManager.getField< fields::acousticvtifields::AcousticLateralSurfaceNodeIndicator >(); + arrayView1d< localIndex > const lateralSurfaceNodeIndicator = nodeManager.getField< acoustivtifields::AcousticLateralSurfaceNodeIndicator >(); /// array of indicators: 1 if a face is on on bottom surface; 0 otherwise - arrayView1d< localIndex > const bottomSurfaceFaceIndicator = faceManager.getField< fields::acousticvtifields::AcousticBottomSurfaceFaceIndicator >(); + arrayView1d< localIndex > const bottomSurfaceFaceIndicator = faceManager.getField< acoustivtifields::AcousticBottomSurfaceFaceIndicator >(); /// array of indicators: 1 if a node is on on bottom surface; 0 otherwise - arrayView1d< localIndex > const bottomSurfaceNodeIndicator = nodeManager.getField< fields::acousticvtifields::AcousticBottomSurfaceNodeIndicator >(); + arrayView1d< localIndex > const bottomSurfaceNodeIndicator = nodeManager.getField< acoustivtifields::AcousticBottomSurfaceNodeIndicator >(); // Lateral surfaces fsManager.apply< FaceManager >( time, @@ -526,9 +490,6 @@ void AcousticVTIZhangWaveEquationSEM::applyFreeSurfaceBC( real64 time, DomainPar /// array of indicators: 1 if a node is on on free surface; 0 otherwise arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfields::AcousticFreeSurfaceNodeIndicator >(); - // freeSurfaceFaceIndicator.zero(); - // freeSurfaceNodeIndicator.zero(); - fsManager.apply< FaceManager >( time, domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), string( "FreeSurface" ), @@ -572,37 +533,120 @@ void AcousticVTIZhangWaveEquationSEM::applyFreeSurfaceBC( real64 time, DomainPar } ); } -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp void AcousticVTIZhangWaveEquationSEM::initializePML() -{ - GEOS_ERROR( "This option is not supported yet" ); +{ GEOS_ERROR( "This option (PML) is not supported" ); return; } -void AcousticVTIZhangWaveEquationSEM::applyPML( real64 const GEOS_UNUSED_PARAM( time ), DomainPartition & GEOS_UNUSED_PARAM( domain ) ) +void AcousticVTIZhangWaveEquationSEM::applyPML( real64 const time, DomainPartition & domain ) { - GEOS_ERROR( "This option is not supported yet" ); - return; + GEOS_MARK_FUNCTION; + + FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); + parametersPML const & param = getReference< parametersPML >( viewKeyStruct::parametersPMLString() ); + + /// Loop over the different mesh bodies; for wave propagation, there is only one mesh body + /// which is the whole mesh + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & ) + { + + NodeManager & nodeManager = mesh.getNodeManager(); + + /// Array views of the pressure p, PML auxiliary variables, and node coordinates + arrayView1d< real32 const > const p_n = nodeManager.getField< acousticfields::Pressure_n >(); + arrayView2d< real32 const > const v_n = nodeManager.getField< acousticfields::AuxiliaryVar1PML >(); + arrayView2d< real32 > const grad_n = nodeManager.getField< acousticfields::AuxiliaryVar2PML >(); + arrayView1d< real32 > const divV_n = nodeManager.getField< acousticfields::AuxiliaryVar3PML >(); + arrayView1d< real32 const > const u_n = nodeManager.getField< acousticfields::AuxiliaryVar4PML >(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords32 = nodeManager.getField< fields::referencePosition32 >().toViewConst(); + + /// Select the subregions concerned by the PML (specified in the xml by the Field Specification) + /// 'targetSet' contains the indices of the elements in a given subregion + fsManager.apply< ElementSubRegionBase, + PerfectlyMatchedLayer >( time, + mesh, + PerfectlyMatchedLayer::catalogName(), + [&]( PerfectlyMatchedLayer const &, + string const &, + SortedArrayView< localIndex const > const & targetSet, + ElementSubRegionBase & subRegion, + string const & ) + + { + + /// Get the element to nodes mapping in the subregion + CellElementSubRegion::NodeMapType const & elemToNodes = + subRegion.getReference< CellElementSubRegion::NodeMapType >( CellElementSubRegion::viewKeyStruct::nodeListString() ); + + /// Get a const ArrayView of the mapping above + traits::ViewTypeConst< CellElementSubRegion::NodeMapType > const elemToNodesViewConst = elemToNodes.toViewConst(); + + /// Array view of the wave speed + arrayView1d< real32 const > const vel = subRegion.getReference< array1d< real32 > >( acousticfields::AcousticVelocity::key()); + + /// Get the object needed to determine the type of the element in the subregion + finiteElement::FiniteElementBase const & + fe = subRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + + real32 xMin[3]; + real32 xMax[3]; + real32 dMin[3]; + real32 dMax[3]; + real32 cMin[3]; + real32 cMax[3]; + for( integer i=0; i<3; ++i ) + { + xMin[i] = param.xMinPML[i]; + xMax[i] = param.xMaxPML[i]; + dMin[i] = param.thicknessMinXYZPML[i]; + dMax[i] = param.thicknessMaxXYZPML[i]; + cMin[i] = param.waveSpeedMinXYZPML[i]; + cMax[i] = param.waveSpeedMaxXYZPML[i]; + } + real32 const r = param.reflectivityPML; + + /// Get the type of the elements in the subregion + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) + { + using FE_TYPE = TYPEOFREF( finiteElement ); + + /// apply the PML kernel + AcousticPMLSEM:: + PMLKernel< FE_TYPE > kernel( finiteElement ); + kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > + ( targetSet, + nodeCoords32, + elemToNodesViewConst, + vel, + p_n, + v_n, + u_n, + xMin, + xMax, + dMin, + dMax, + cMin, + cMax, + r, + grad_n, + divV_n ); + } ); + } ); + } ); + } real64 AcousticVTIZhangWaveEquationSEM::explicitStepForward( real64 const & time_n, - real64 const & dt, - integer cycleNumber, - DomainPartition & domain, - bool computeGradient ) -{ - real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain, true ); -======= -real64 AcousticVTIWaveEquationSEM::explicitStepForward( real64 const & time_n, - real64 const & dt, - integer, - DomainPartition & domain, - bool computeGradient ) + real64 const & dt, + integer cycleNumber, + DomainPartition & domain, + bool computeGradient ) { - real64 dtOut = explicitStepInternal( time_n, dt, domain ); ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp + real64 dtCompute = explicitStepInternal( time_n, dt, domain, true ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -627,40 +671,26 @@ real64 AcousticVTIWaveEquationSEM::explicitStepForward( real64 const & time_n, prepareNextTimestep( mesh ); } ); - return dtOut; + return dtCompute; } - real64 AcousticVTIZhangWaveEquationSEM::explicitStepBackward( real64 const & time_n, - real64 const & dt, - integer cycleNumber, - DomainPartition & domain, - bool computeGradient ) -{ -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp - real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain, false ); -======= - GEOS_ERROR( "This option is not supported yet" ); - return -1; -} - -real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, - real64 const & dt, - DomainPartition & domain ) + real64 const & dt, + integer cycleNumber, + DomainPartition & domain, + bool computeGradient ) { - GEOS_MARK_FUNCTION; - - GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc ); ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp - + real64 dtCompute = explicitStepInternal( time_n, dt, domain, false ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, - arrayView1d< string const > const & GEOS_UNUSED_PARAM ( regionNames ) ) + arrayView1d< string const > const & regionNames ) { NodeManager & nodeManager = mesh.getNodeManager(); + arrayView1d< real32 const > const mass = nodeManager.getField< acousticfields::AcousticMassVector >(); + arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticvtifields::Pressure_p_nm1 >(); arrayView1d< real32 > const p_n = nodeManager.getField< acousticvtifields::Pressure_p_n >(); arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticvtifields::Pressure_p_np1 >(); @@ -669,6 +699,11 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, arrayView1d< real32 > const q_n = nodeManager.getField< acousticvtifields::Pressure_q_n >(); arrayView1d< real32 > const q_np1 = nodeManager.getField< acousticvtifields::Pressure_q_np1 >(); + + EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); + real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() ); + int const maxCycle = int(round( maxTime / dt )); + if( computeGradient && cycleNumber >= 0 ) { GEOS_ERROR( "This option is not supported yet" ); @@ -677,7 +712,7 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, prepareNextTimestep( mesh ); } ); - return dtOut; + return dtCompute; } void AcousticVTIZhangWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh ) @@ -713,13 +748,12 @@ void AcousticVTIZhangWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh ) } ); } -void AcousticVTIZhangWaveEquationSEM::computeUnknowns( real64 const & GEOS_UNUSED_PARAM( time_n ), - real64 const & dt, - integer cycleNumber, - DomainPartition & GEOS_UNUSED_PARAM( domain ), - MeshLevel & mesh, - arrayView1d< string const > const & regionNames, - bool const isForward ) +void AcousticVTIZhangWaveEquationSEM::computeUnknowns( real64 const & time_n, + real64 const & dt, + DomainPartition & domain, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames, + bool const isForward ) { NodeManager & nodeManager = mesh.getNodeManager(); @@ -737,14 +771,17 @@ void AcousticVTIZhangWaveEquationSEM::computeUnknowns( real64 const & GEOS_UNUSE arrayView1d< real32 > const q_n = nodeManager.getField< acousticvtifields::Pressure_q_n >(); arrayView1d< real32 > const q_np1 = nodeManager.getField< acousticvtifields::Pressure_q_np1 >(); + arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticfields::Pressure_nm1 >(); + arrayView1d< real32 > const p_n = nodeManager.getField< acousticfields::Pressure_n >(); + arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticfields::Pressure_np1 >(); + arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfields::AcousticFreeSurfaceNodeIndicator >(); arrayView1d< localIndex const > const lateralSurfaceNodeIndicator = nodeManager.getField< acousticvtifields::AcousticLateralSurfaceNodeIndicator >(); arrayView1d< localIndex const > const bottomSurfaceNodeIndicator = nodeManager.getField< acousticvtifields::AcousticBottomSurfaceNodeIndicator >(); arrayView1d< real32 > const stiffnessVector_p = nodeManager.getField< acousticvtifields::StiffnessVector_p >(); arrayView1d< real32 > const stiffnessVector_q = nodeManager.getField< acousticvtifields::StiffnessVector_q >(); arrayView1d< real32 > const rhs = nodeManager.getField< acousticfields::ForcingRHS >(); - - if( isForward ) +if( isForward ) { auto kernelFactory = acousticVTIZhangWaveEquationSEMKernels::ExplicitAcousticVTIZhangSEMFactory( dt ); @@ -770,27 +807,19 @@ void AcousticVTIZhangWaveEquationSEM::computeUnknowns( real64 const & GEOS_UNUSE "", kernelFactory ); -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp } //Modification of cycleNember useful when minTime < 0 - EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); - real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() ); - integer const cycleForSource = int(round( -minTime / dt + cycleNumber )); - addSourceToRightHandSide( cycleForSource, rhs ); -======= - addSourceToRightHandSide( time_n, rhs ); + addSourceToRightHandSide( time_n, rhs ); - /// calculate your time integrators ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp + /// calculate your time integrators + real64 const dt2 = pow( dt, 2 ); SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst(); if( !m_usePML ) { GEOS_MARK_SCOPE ( updateP ); - AcousticTimeSchemeSEM::LeapFrogVTIWithoutPML( dt, p_np1, p_n, p_nm1, q_np1, q_n, q_nm1, mass, stiffnessVector_p, - stiffnessVector_q, damping_pp, damping_pq, damping_qp, damping_qq, - rhs, freeSurfaceNodeIndicator, lateralSurfaceNodeIndicator, - bottomSurfaceNodeIndicator, solverTargetNodesSet ); + AcousticTimeSchemeSEM::LeapFrogWithoutPML( dt, p_np1, p_n, p_nm1, mass, stiffnessVector, damping, + rhs, freeSurfaceNodeIndicator, solverTargetNodesSet ); } else { @@ -799,11 +828,10 @@ void AcousticVTIZhangWaveEquationSEM::computeUnknowns( real64 const & GEOS_UNUSE } void AcousticVTIZhangWaveEquationSEM::synchronizeUnknowns( real64 const & time_n, - real64 const & dt, - integer const, - DomainPartition & domain, - MeshLevel & mesh, - arrayView1d< string const > const & ) + real64 const & dt, + DomainPartition & domain, + MeshLevel & mesh, + arrayView1d< string const > const & ) { NodeManager & nodeManager = mesh.getNodeManager(); @@ -826,6 +854,7 @@ void AcousticVTIZhangWaveEquationSEM::synchronizeUnknowns( real64 const & time_n GEOS_ERROR( "This option is not supported yet" ); } + CommunicationTools & syncFields = CommunicationTools::getInstance(); syncFields.synchronizeFields( fieldsToBeSync, mesh, @@ -845,31 +874,33 @@ void AcousticVTIZhangWaveEquationSEM::synchronizeUnknowns( real64 const & time_n } real64 AcousticVTIZhangWaveEquationSEM::explicitStepInternal( real64 const & time_n, - real64 const & dt, - integer const cycleNumber, - DomainPartition & domain, - bool const isForward ) + real64 const & dt, + DomainPartition & domain, + bool const isForward ) { GEOS_MARK_FUNCTION; GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc ); + real64 dtCompute; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - computeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames, isForward ); - synchronizeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames ); + localIndex nSubSteps = (int) ceil( dt/m_timeStep ); + dtCompute = dt/nSubSteps; + computeUnknowns( time_n, dtCompute, domain, mesh, regionNames, isForward ); + synchronizeUnknowns( time_n, dtCompute, domain, mesh, regionNames ); } ); - return dt; + return dtCompute; } void AcousticVTIZhangWaveEquationSEM::cleanup( real64 const time_n, - integer const cycleNumber, - integer const eventCounter, - real64 const eventProgress, - DomainPartition & domain ) + integer const cycleNumber, + integer const eventCounter, + real64 const eventProgress, + DomainPartition & domain ) { // call the base class cleanup (for reporting purposes) PhysicsSolverBase::cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain ); @@ -893,10 +924,6 @@ void AcousticVTIZhangWaveEquationSEM::cleanup( real64 const time_n, } ); } -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp -REGISTER_CATALOG_ENTRY( SolverBase, AcousticVTIZhangWaveEquationSEM, string const &, dataRepository::Group * const ) -======= -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, AcousticVTIWaveEquationSEM, string const &, dataRepository::Group * const ) ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, AcousticVTIZhangWaveEquationSEM, string const &, dataRepository::Group * const ) } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.hpp index 0096d80240f..d9dc6a9ffb4 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.hpp @@ -23,12 +23,7 @@ #include "physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp" #include "mesh/MeshFields.hpp" -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.hpp -#include "physicsSolvers/SolverBase.hpp" -======= #include "physicsSolvers/PhysicsSolverBase.hpp" -#include "physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp" ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp" #include "AcousticVTIFields.hpp" @@ -43,7 +38,7 @@ class AcousticVTIZhangWaveEquationSEM : public WaveSolverBase using ATOMIC_POLICY = AtomicPolicy< EXEC_POLICY >; AcousticVTIZhangWaveEquationSEM( const std::string & name, - Group * const parent ); + Group * const parent ); virtual ~AcousticVTIZhangWaveEquationSEM() override; @@ -101,6 +96,11 @@ class AcousticVTIZhangWaveEquationSEM : public WaveSolverBase */ virtual void initializePML() override; + /** + */ + virtual real64 computeTimeStep( real64 & dtOut ) override; + + /** * @brief Overridden from ExecutableGroup. Used to write last seismogram if needed. @@ -124,17 +124,11 @@ class AcousticVTIZhangWaveEquationSEM : public WaveSolverBase */ real64 explicitStepInternal( real64 const & time_n, real64 const & dt, -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.hpp - integer const cycleNumber, DomainPartition & domain, bool const isForward ); -======= - DomainPartition & domain ); ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp void computeUnknowns( real64 const & time_n, real64 const & dt, - integer const cycleNumber, DomainPartition & domain, MeshLevel & mesh, arrayView1d< string const > const & regionNames, @@ -142,7 +136,6 @@ class AcousticVTIZhangWaveEquationSEM : public WaveSolverBase void synchronizeUnknowns( real64 const & time_n, real64 const & dt, - integer const cycleNumber, DomainPartition & domain, MeshLevel & mesh, arrayView1d< string const > const & regionNames ); @@ -165,12 +158,6 @@ class AcousticVTIZhangWaveEquationSEM : public WaveSolverBase */ virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - /** - * @brief Compute the lateral and bottom surface Field indicators of the boxed domain - * @param domain the partition domain - */ - virtual void precomputeSurfaceFieldIndicator( DomainPartition & domain ); - /** * @brief Apply free surface condition to the face define in the geometry box from the xml * @param time the time to apply the BC @@ -178,7 +165,6 @@ class AcousticVTIZhangWaveEquationSEM : public WaveSolverBase */ virtual void applyFreeSurfaceBC( real64 const time, DomainPartition & domain ) override; -<<<<<<< HEAD:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.hpp /** * @brief Apply Perfectly Matched Layer (PML) to the regions defined in the geometry box from the xml * @param time the time to apply the BC @@ -187,14 +173,9 @@ class AcousticVTIZhangWaveEquationSEM : public WaveSolverBase virtual void applyPML( real64 const time, DomainPartition & domain ) override; /// Pressure_np1 at the receiver location for each time step for each receiver -======= - virtual real64 computeTimeStep( real64 & dtOut ) override; - - /// Pressure_p_np1 at the receiver location for each time step for each receiver ->>>>>>> develop:src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp array2d< real32 > m_pressureNp1AtReceivers; - /// Array of size the number of receivers and full of 0.5 (used for calculating the seismos) + /// Array of size the number of receivers and filled with 0.5 (used for calculating the seismos) array1d< real32 > m_seismoCoeff; }; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEMKernel.hpp index 9e128110687..c81d0659c1a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEMKernel.hpp @@ -26,18 +26,17 @@ #include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" #endif #include "physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp" +#include "AcousticVTIFields.hpp" namespace geos { -using namespace fields; - /// Namespace to contain the acoustic wave kernels. namespace acousticVTIZhangWaveEquationSEMKernels { /** - * @brief Implements kernels for solving the acoustic wave equations + * @brief Implements kernels for solving the pseudo-acoustic VTI wave equations * explicit central FD method and SEM * @copydoc geos::finiteElement::KernelBase * @tparam SUBREGION_TYPE The type of subregion that the kernel will act on. @@ -56,11 +55,11 @@ namespace acousticVTIZhangWaveEquationSEMKernels template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > -class ExplicitAcousticVTIZhangSEM : public finiteElement::KernelBase< SUBREGION_TYPE, - CONSTITUTIVE_TYPE, - FE_TYPE, - 1, - 1 > +class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 1, + 1 > { public: @@ -94,14 +93,14 @@ class ExplicitAcousticVTIZhangSEM : public finiteElement::KernelBase< SUBREGION_ * @param dt The time interval for the step. * elements to be processed during this kernel launch. */ - ExplicitAcousticVTIZhangSEM( NodeManager & nodeManager, - EdgeManager const & edgeManager, - FaceManager const & faceManager, - localIndex const targetRegionIndex, - SUBREGION_TYPE const & elementSubRegion, - FE_TYPE const & finiteElementSpace, - CONSTITUTIVE_TYPE & inputConstitutiveType, - real64 const dt ): + ExplicitAcousticSEM( NodeManager & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + real64 const dt ): Base( elementSubRegion, finiteElementSpace, inputConstitutiveType ), @@ -124,7 +123,7 @@ class ExplicitAcousticVTIZhangSEM : public finiteElement::KernelBase< SUBREGION_ /** * @copydoc geos::finiteElement::KernelBase::StackVariables * - * ### ExplicitAcousticVTIZhangSEM Description + * ### ExplicitAcousticSEM Description * Adds a stack arrays for the nodal force, primary displacement variable, etc. */ struct StackVariables : Base::StackVariables @@ -199,7 +198,7 @@ class ExplicitAcousticVTIZhangSEM : public finiteElement::KernelBase< SUBREGION_ /** * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel * - * ### ExplicitAcousticVTIZhangSEM Description + * ### ExplicitAcousticSEM Description * Calculates stiffness vector * */ @@ -209,21 +208,22 @@ class ExplicitAcousticVTIZhangSEM : public finiteElement::KernelBase< SUBREGION_ localIndex const q, StackVariables & stack ) const { - // Pseudo Stiffness xy + // (A_xy nabla u)((A_xy nabla v) m_finiteElementSpace.template computeStiffnessxyTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_p = -val * stack.invDensity * stack.vti_epsi * m_p_n[m_elemsToNodes( k, j )]; stack.stiffnessVectorLocal_p[i] += localIncrement_p; + real32 const localIncrement_q = -val * stack.invDensity * stack.vti_sqrtDelta * m_p_n[m_elemsToNodes( k, j )]; stack.stiffnessVectorLocal_q[i] += localIncrement_q; } ); - // Pseudo-Stiffness z - + // (A_z nabla u)((A_z nabla v) m_finiteElementSpace.template computeStiffnesszTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) { real32 const localIncrement_p = -val * stack.invDensity * stack.vti_sqrtDelta* m_q_n[m_elemsToNodes( k, j )]; stack.stiffnessVectorLocal_p[i] += localIncrement_p; + real32 const localIncrement_q = -val * stack.invDensity * m_q_n[m_elemsToNodes( k, j )]; stack.stiffnessVectorLocal_q[i] += localIncrement_q; } ); @@ -260,9 +260,9 @@ class ExplicitAcousticVTIZhangSEM : public finiteElement::KernelBase< SUBREGION_ -/// The factory used to construct a ExplicitAcousticWaveEquation kernel. +/// The factory used to construct a ExplicitAcousticVTIZhang kernel. using ExplicitAcousticVTIZhangSEMFactory = finiteElement::KernelFactory< ExplicitAcousticVTIZhangSEM, - real64 >; + real64 >; } // namespace acousticVTIZhangWaveEquationSEMKernels