diff --git a/integratedTests b/integratedTests index c86fa430fb1..a9591c65638 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit c86fa430fb1c676ded5b10ca3294dd405ad3f73b +Subproject commit a9591c656385d5af804e6292c0f4ab0dfd770f42 diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp index 90e66c66167..c7ec3d9d330 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp @@ -65,21 +65,10 @@ AcousticFirstOrderWaveEquationSEM::AcousticFirstOrderWaveEquationSEM( const std: setSizedFromParent( 0 ). setDescription( "Element containing the sources" ); - registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Element containing the receivers" ); - registerWrapper( viewKeyStruct::sourceRegionString(), &m_sourceRegion ). setInputFlag( InputFlags::FALSE ). setSizedFromParent( 0 ). setDescription( "Region containing the sources" ); - - registerWrapper( viewKeyStruct::receiverRegionString(), &m_receiverRegion ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Region containing the receivers" ); - } AcousticFirstOrderWaveEquationSEM::~AcousticFirstOrderWaveEquationSEM() @@ -151,10 +140,10 @@ void AcousticFirstOrderWaveEquationSEM::postProcessInput() localIndex const numSourcesGlobal = m_sourceCoordinates.size( 0 ); localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 ); - m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_uxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_uyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_uzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); + m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_uxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_uyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_uzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); m_sourceElem.resize( numSourcesGlobal ); m_sourceRegion.resize( numSourcesGlobal ); m_rcvElem.resize( numReceiversGlobal ); @@ -288,8 +277,7 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); - real64 const time = 0.0; - applyFreeSurfaceBC( time, domain ); + applyFreeSurfaceBC( 0.0, domain ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, @@ -354,6 +342,8 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro } ); } ); } ); + + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); } @@ -524,9 +514,9 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t cycleNumber, p_np1 ); } ); - arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView(); - arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView(); - arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView(); + arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView(); + arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView(); + arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView(); compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_x, velocity_x, uxReceivers ); compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_y, velocity_y, uyReceivers ); @@ -545,16 +535,12 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t true ); // compute the seismic traces since last step. - arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); computeAllSeismoTraces( time_n, dt, p_np1, p_np1, pReceivers ); - // increment m_indexSeismoTrace - while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) - { - m_indexSeismoTrace++; - } - + incrementIndexSeismoTrace( time_n ); } ); + return dt; } @@ -574,62 +560,26 @@ void AcousticFirstOrderWaveEquationSEM::cleanup( real64 const time_n, integer co arrayView2d< real32 > const velocity_y = elementSubRegion.getField< wavesolverfields::Velocity_y >(); arrayView2d< real32 > const velocity_z = elementSubRegion.getField< wavesolverfields::Velocity_z >(); - arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView(); - arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView(); - arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView(); - - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_x, velocity_x, uxReceivers ); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_y, velocity_y, uyReceivers ); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_z, velocity_z, uzReceivers ); + arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView(); + arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView(); + arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView(); - } ); - arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); - computeAllSeismoTraces( time_n, 0, p_np1, p_np1, pReceivers ); - } ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_x, velocity_x, uxReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_y, velocity_y, uyReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_z, velocity_z, uzReceivers ); - // increment m_indexSeismoTrace - while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) - { - m_indexSeismoTrace++; - } -} + WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, uxReceivers, uyReceivers, uzReceivers ); -void AcousticFirstOrderWaveEquationSEM::computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ) -{ - localIndex indexSeismoTrace = m_indexSeismoTrace; - for( real64 timeSeismo; - (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace; - indexSeismoTrace++ ) - { - WaveSolverUtils::computeSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal, m_nsamplesSeismoTrace, - m_outputSeismoTrace, - var_np1, var_n, varAtReceivers ); - } -} + } ); + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); + computeAllSeismoTraces( time_n, 0.0, p_np1, p_np1, pReceivers ); + WaveSolverUtils::writeSeismoTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, pReceivers ); -void AcousticFirstOrderWaveEquationSEM::compute2dVariableAllSeismoTraces( localIndex const regionIndex, - real64 const time_n, - real64 const dt, - arrayView2d< real32 const > const var_np1, - arrayView2d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ) -{ - localIndex indexSeismoTrace = m_indexSeismoTrace; - for( real64 timeSeismo; - (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace; - indexSeismoTrace++ ) - { - WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dt, regionIndex, m_receiverRegion, timeSeismo, indexSeismoTrace, m_rcvElem, m_receiverConstants, m_receiverIsLocal, m_nsamplesSeismoTrace, - m_outputSeismoTrace, - var_np1, var_n, varAtReceivers ); - } + } ); } - void AcousticFirstOrderWaveEquationSEM::initializePML() { GEOS_ERROR( getDataContext() << ": PML for the first order acoustic wave propagator not yet implemented" ); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp index 97bcfe470ce..709d5f86789 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp @@ -21,8 +21,8 @@ #define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICFIRSTORDERWAVEEQUATIONSEM_HPP_ #include "mesh/MeshFields.hpp" -#include "WaveSolverUtils.hpp" #include "WaveSolverBaseFields.hpp" +#include "WaveSolverBase.hpp" namespace geos { @@ -78,37 +78,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase */ virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ); - /** - * TODO: move implementation into WaveSolverUtils - * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt - * @param time_n the time corresponding to the field values pressure_n - * @param dt the simulation timestep - * @param var_np1 the field values at time_n + dt - * @param var_n the field values at time_n - * @param var_receivers the array holding the trace values, where the output is written - */ - virtual void computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ); - - /** - * TODO: move implementation into WaveSolverUtils - * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt for a 2d variable - * @param time_n the time corresponding to the field values pressure_n - * @param dt the simulation timestep - * @param var_np1 the field values at time_n + dt - * @param var_n the field values at time_n - * @param var_receivers the array holding the trace values, where the output is written - */ - virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex, - real64 const time_n, - real64 const dt, - arrayView2d< real32 const > const var_np1, - arrayView2d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ); - /** * @brief Initialize Perfectly Matched Layer (PML) information @@ -133,8 +102,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase static constexpr char const * sourceElemString() { return "sourceElem"; } static constexpr char const * sourceRegionString() { return "sourceRegion"; } - static constexpr char const * receiverElemString() { return "rcvElem"; } - static constexpr char const * receiverRegionString() { return "receiverRegion"; } } waveEquationViewKeys; @@ -197,13 +164,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase /// Array containing the elements which contain the region which the source belongs array1d< localIndex > m_sourceRegion; - - /// Array containing the elements which contain the region which the receiver belongs - array1d< localIndex > m_receiverRegion; - - /// Array containing the elements which contain a receiver - array1d< localIndex > m_rcvElem; - }; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp index 68dff0e5793..8922a92e8cd 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp @@ -106,7 +106,7 @@ void AcousticVTIWaveEquationSEM::postProcessInput() localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 ); - m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); + m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); } void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, @@ -229,8 +229,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups() DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); - real64 const time = 0.0; - applyFreeSurfaceBC( time, domain ); + applyFreeSurfaceBC( 0.0, domain ); precomputeSurfaceFieldIndicator( domain ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -315,6 +314,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); } void AcousticVTIWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain ) @@ -649,10 +649,11 @@ real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, true ); // compute the seismic traces since last step. - arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); - + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers ); + incrementIndexSeismoTrace( time_n ); + /// prepare next step forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) { @@ -681,45 +682,16 @@ void AcousticVTIWaveEquationSEM::cleanup( real64 const time_n, arrayView1d< string const > const & ) { NodeManager & nodeManager = mesh.getNodeManager(); - arrayView1d< real32 const > const p_n = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >(); + arrayView1d< real32 const > const p_n = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >(); arrayView1d< real32 const > const p_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_np1 >(); - arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); - computeAllSeismoTraces( time_n, 0, p_np1, p_n, pReceivers ); - } ); -} + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); + computeAllSeismoTraces( time_n, 0.0, p_np1, p_n, pReceivers ); -void AcousticVTIWaveEquationSEM::computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ) -{ - - /* - * In forward case we compute seismo if time_n + dt is the first time - * step after the timeSeismo to write. - * - * time_n timeSeismo time_n + dt - * ---|--------------|-------------| - * - * In backward (time_n goes decreasing) case we compute seismo if - * time_n is the last time step before the timeSeismo to write. - * - * time_n - dt timeSeismo time_n - * ---|--------------|-------------| - */ - for( real64 timeSeismo; - (m_forward)?((timeSeismo = m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + dt + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace): - ((timeSeismo = m_dtSeismoTrace*(m_nsamplesSeismoTrace-m_indexSeismoTrace-1)) >= (time_n - dt - epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace); - m_indexSeismoTrace++ ) - { - WaveSolverUtils::computeSeismoTrace( time_n, (m_forward)?dt:-dt, timeSeismo, (m_forward)?m_indexSeismoTrace:(m_nsamplesSeismoTrace-m_indexSeismoTrace-1), m_receiverNodeIds, m_receiverConstants, - m_receiverIsLocal, - m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers ); - } + WaveSolverUtils::writeSeismoTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, pReceivers ); + } ); } - REGISTER_CATALOG_ENTRY( SolverBase, AcousticVTIWaveEquationSEM, string const &, dataRepository::Group * const ) } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp index 1fe0ede62af..258cbbebb64 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp @@ -72,21 +72,6 @@ class AcousticVTIWaveEquationSEM : public WaveSolverBase */ virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ); - /** - * TODO: move implementation into WaveSolverBase - * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt - * @param time_n the time corresponding to the field values pressure_n - * @param dt the simulation timestep - * @param var_np1 the field values at time_n + dt - * @param var_n the field values at time_n - * @param varAtReceivers the array holding the trace values, where the output is written - */ - virtual void computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ); - /** * @brief Overridden from ExecutableGroup. Used to write last seismogram if needed. */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp index 04e18f4701b..6cbadc9ec14 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp @@ -109,10 +109,7 @@ void AcousticWaveEquationSEM::postProcessInput() WaveSolverBase::postProcessInput(); - localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 ); - - m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - + m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, m_receiverCoordinates.size( 0 ) + 1 ); } void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, @@ -248,8 +245,7 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); - real64 const time = 0.0; - applyFreeSurfaceBC( time, domain ); + applyFreeSurfaceBC( 0.0, domain ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, @@ -327,6 +323,7 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); } @@ -1067,11 +1064,10 @@ real64 AcousticWaveEquationSEM::explicitStepInternal( real64 const & time_n, true ); /// compute the seismic traces since last step. - arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); - if( time_n >= 0 ) - { - computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers ); - } + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); + computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers ); + incrementIndexSeismoTrace( time_n ); + /// prepare next step forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) { @@ -1109,40 +1105,12 @@ void AcousticWaveEquationSEM::cleanup( real64 const time_n, NodeManager & nodeManager = mesh.getNodeManager(); arrayView1d< real32 const > const p_n = nodeManager.getField< fields::Pressure_n >(); arrayView1d< real32 const > const p_np1 = nodeManager.getField< fields::Pressure_np1 >(); - arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); - computeAllSeismoTraces( time_n, 0, p_np1, p_n, pReceivers ); - } ); -} - -void AcousticWaveEquationSEM::computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ) -{ + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); + computeAllSeismoTraces( time_n, 0.0, p_np1, p_n, pReceivers ); - /* - * In forward case we compute seismo if time_n + dt is the first time - * step after the timeSeismo to write. - * - * time_n timeSeismo time_n + dt - * ---|--------------|-------------| - * - * In backward (time_n goes decreasing) case we compute seismo if - * time_n is the last time step before the timeSeismo to write. - * - * time_n - dt timeSeismo time_n - * ---|--------------|-------------| - */ - for( real64 timeSeismo; - (m_forward)?((timeSeismo = m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + dt + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace): - ((timeSeismo = m_dtSeismoTrace*(m_nsamplesSeismoTrace-m_indexSeismoTrace-1)) >= (time_n - dt - epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace); - m_indexSeismoTrace++ ) - { - WaveSolverUtils::computeSeismoTrace( time_n, (m_forward)?dt:-dt, timeSeismo, (m_forward)?m_indexSeismoTrace:(m_nsamplesSeismoTrace-m_indexSeismoTrace-1), m_receiverNodeIds, m_receiverConstants, - m_receiverIsLocal, - m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers ); - } + WaveSolverUtils::writeSeismoTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, pReceivers ); + } ); } REGISTER_CATALOG_ENTRY( SolverBase, AcousticWaveEquationSEM, string const &, dataRepository::Group * const ) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp index b30bd7cf6c6..fa925902693 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp @@ -81,21 +81,6 @@ class AcousticWaveEquationSEM : public WaveSolverBase */ virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ); - /** - * TODO: move implementation into WaveSolverBase - * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt - * @param time_n the time corresponding to the field values pressure_n - * @param dt the simulation timestep - * @param var_np1 the field values at time_n + dt - * @param var_n the field values at time_n - * @param varAtReceivers the array holding the trace values, where the output is written - */ - virtual void computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ); - /** * @brief Initialize Perfectly Matched Layer (PML) information @@ -110,14 +95,6 @@ class AcousticWaveEquationSEM : public WaveSolverBase struct viewKeyStruct : WaveSolverBase::viewKeyStruct { - static constexpr char const * sourceNodeIdsString() { return "sourceNodeIds"; } - static constexpr char const * sourceConstantsString() { return "sourceConstants"; } - static constexpr char const * sourceIsAccessibleString() { return "sourceIsAccessible"; } - - static constexpr char const * receiverNodeIdsString() { return "receiverNodeIds"; } - static constexpr char const * receiverConstantsString() {return "receiverConstants"; } - static constexpr char const * receiverIsLocalString() { return "receiverIsLocal"; } - static constexpr char const * pressureNp1AtReceiversString() { return "pressureNp1AtReceivers"; } } waveEquationViewKeys; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp index e0d673af714..1e4364f840c 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp @@ -89,21 +89,10 @@ ElasticFirstOrderWaveEquationSEM::ElasticFirstOrderWaveEquationSEM( const std::s setSizedFromParent( 0 ). setDescription( "Element containing the sources" ); - registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Element containing the receivers" ); - registerWrapper( viewKeyStruct::sourceRegionString(), &m_sourceRegion ). setInputFlag( InputFlags::FALSE ). setSizedFromParent( 0 ). setDescription( "Region containing the sources" ); - - registerWrapper( viewKeyStruct::receiverRegionString(), &m_receiverRegion ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Region containing the receivers" ); - } ElasticFirstOrderWaveEquationSEM::~ElasticFirstOrderWaveEquationSEM() @@ -196,18 +185,16 @@ void ElasticFirstOrderWaveEquationSEM::postProcessInput() m_rcvElem.resize( numReceiversGlobal ); m_receiverRegion.resize( numReceiversGlobal ); - m_displacementxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_displacementyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_displacementzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - - m_sigmaxxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_sigmayyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_sigmazzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_sigmaxyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_sigmaxzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_sigmayzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - + m_displacementxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_displacementyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_displacementzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_sigmaxxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_sigmayyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_sigmazzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_sigmaxyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_sigmaxzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_sigmayzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); } @@ -642,8 +629,6 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressxy, stressxy, sigmaxyReceivers ); compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressxz, stressxz, sigmaxzReceivers ); compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressyz, stressyz, sigmayzReceivers ); - - } ); FieldIdentifiers fieldsToBeSync; @@ -668,12 +653,7 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti computeAllSeismoTraces( time_n, dt, uy_np1, uy_np1, uyReceivers ); computeAllSeismoTraces( time_n, dt, uz_np1, uz_np1, uzReceivers ); - // increment m_indexSeismoTrace - while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) - { - m_indexSeismoTrace++; - } - + incrementIndexSeismoTrace( time_n ); } ); return dt; @@ -712,68 +692,35 @@ void ElasticFirstOrderWaveEquationSEM::cleanup( real64 const time_n, arrayView2d< real32 > const sigmaxzReceivers = m_sigmaxzNp1AtReceivers.toView(); arrayView2d< real32 > const sigmayzReceivers = m_sigmayzNp1AtReceivers.toView(); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressxx, stressxx, sigmaxxReceivers ); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressyy, stressyy, sigmayyReceivers ); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stresszz, stresszz, sigmazzReceivers ); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressxy, stressxy, sigmaxyReceivers ); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressxz, stressxz, sigmaxzReceivers ); - compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressyz, stressyz, sigmayzReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressxx, stressxx, sigmaxxReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressyy, stressyy, sigmayyReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stresszz, stresszz, sigmazzReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressxy, stressxy, sigmaxyReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressxz, stressxz, sigmaxzReceivers ); + compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressyz, stressyz, sigmayzReceivers ); + + WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, sigmaxxReceivers, sigmayyReceivers, sigmazzReceivers ); + WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, sigmaxyReceivers, sigmaxzReceivers, sigmayzReceivers ); + } ); arrayView1d< real32 > const ux_np1 = nodeManager.getField< wavesolverfields::Displacementx_np1 >(); arrayView1d< real32 > const uy_np1 = nodeManager.getField< wavesolverfields::Displacementy_np1 >(); arrayView1d< real32 > const uz_np1 = nodeManager.getField< wavesolverfields::Displacementz_np1 >(); // compute the seismic traces since last step. - arrayView2d< real32 > const uxReceivers = m_displacementxNp1AtReceivers.toView(); - arrayView2d< real32 > const uyReceivers = m_displacementyNp1AtReceivers.toView(); - arrayView2d< real32 > const uzReceivers = m_displacementzNp1AtReceivers.toView(); + arrayView2d< real32 > const uxReceivers = m_displacementxNp1AtReceivers.toView(); + arrayView2d< real32 > const uyReceivers = m_displacementyNp1AtReceivers.toView(); + arrayView2d< real32 > const uzReceivers = m_displacementzNp1AtReceivers.toView(); - computeAllSeismoTraces( time_n, 0, ux_np1, ux_np1, uxReceivers ); - computeAllSeismoTraces( time_n, 0, uy_np1, uy_np1, uyReceivers ); - computeAllSeismoTraces( time_n, 0, uz_np1, uz_np1, uzReceivers ); + computeAllSeismoTraces( time_n, 0.0, ux_np1, ux_np1, uxReceivers ); + computeAllSeismoTraces( time_n, 0.0, uy_np1, uy_np1, uyReceivers ); + computeAllSeismoTraces( time_n, 0.0, uz_np1, uz_np1, uzReceivers ); + WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, uxReceivers, uyReceivers, uzReceivers ); } ); - -// increment m_indexSeismoTrace - while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) - { - m_indexSeismoTrace++; - } - -} - -void ElasticFirstOrderWaveEquationSEM::computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ) -{ - - localIndex indexSeismoTrace = m_indexSeismoTrace; - for( real64 timeSeismo; - (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace; - indexSeismoTrace++ ) - { - WaveSolverUtils::computeSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal, - m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers ); - } -} - -void ElasticFirstOrderWaveEquationSEM::compute2dVariableAllSeismoTraces( localIndex const regionIndex, - real64 const time_n, - real64 const dt, - arrayView2d< real32 const > const var_np1, - arrayView2d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ) -{ - localIndex indexSeismoTrace = m_indexSeismoTrace; - for( real64 timeSeismo; - (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace; - indexSeismoTrace++ ) - { - WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dt, regionIndex, m_receiverRegion, timeSeismo, indexSeismoTrace, m_rcvElem, m_receiverConstants, m_receiverIsLocal, - m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers ); - } } void ElasticFirstOrderWaveEquationSEM::initializePML() diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp index a4c3c6a3266..c93d7cb4eeb 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp @@ -21,9 +21,8 @@ #define SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ELASTICFIRSTORDERWAVEEQUATIONSEM_HPP_ #include "mesh/MeshFields.hpp" -#include "WaveSolverUtils.hpp" #include "WaveSolverBaseFields.hpp" - +#include "WaveSolverBase.hpp" namespace geos @@ -77,36 +76,6 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase */ void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ); - /** - * TODO: move implementation into WaveSolverBase - * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt - * @param time_n the time corresponding to the field values pressure_n - * @param dt the simulation timestep - * @param var_np1 the field values at time_n + dt - * @param var_n the field values at time_n - * @param varAtreceivers the array holding the trace values, where the output is written - */ - virtual void computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ); - - /** - * TODO: move implementation into WaveSolverBase - * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt - * @param time_n the time corresponding to the field values pressure_n - * @param dt the simulation timestep - * @param var_np1 the field values at time_n + dt - * @param var_n the field values at time_n - * @param varAtreceivers the array holding the trace values, where the output is written - */ - virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex, - real64 const time_n, - real64 const dt, - arrayView2d< real32 const > const var_np1, - arrayView2d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ); /** * @brief Initialize Perfectly Matched Layer (PML) information @@ -135,8 +104,6 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase static constexpr char const * sourceElemString() { return "sourceElem"; } static constexpr char const * sourceRegionString() { return "sourceRegion"; } - static constexpr char const * receiverElemString() { return "rcvElem"; } - static constexpr char const * receiverRegionString() { return "receiverRegion"; } } waveEquationViewKeys; @@ -216,12 +183,6 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase /// Array containing the elements which contain the region which the source belongs array1d< localIndex > m_sourceRegion; - /// Array containing the elements which contain a receiver - array1d< localIndex > m_rcvElem; - - /// Array containing the elements which contain the region which the receiver belongs - array1d< localIndex > m_receiverRegion; - }; } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp index a6f5a214281..2de2021b306 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp @@ -38,12 +38,6 @@ ElasticWaveEquationSEM::ElasticWaveEquationSEM( const std::string & name, WaveSolverBase( name, parent ) { - - registerWrapper( viewKeyStruct::sourceNodeIdsString(), &m_sourceNodeIds ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Indices of the nodes (in the right order) for each source point" ); - registerWrapper( viewKeyStruct::sourceConstantsString(), &m_sourceConstantsx ). setInputFlag( InputFlags::FALSE ). setSizedFromParent( 0 ). @@ -59,26 +53,6 @@ ElasticWaveEquationSEM::ElasticWaveEquationSEM( const std::string & name, setSizedFromParent( 0 ). setDescription( "Constant part of the source for the nodes listed in m_sourceNodeIds in z-direction" ); - registerWrapper( viewKeyStruct::sourceIsAccessibleString(), &m_sourceIsAccessible ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Flag that indicates whether the source is accessible to this MPI rank" ); - - registerWrapper( viewKeyStruct::receiverNodeIdsString(), &m_receiverNodeIds ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Indices of the nodes (in the right order) for each receiver point" ); - - registerWrapper( viewKeyStruct::sourceConstantsString(), &m_receiverConstants ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Constant part of the receiver for the nodes listed in m_receiverNodeIds" ); - - registerWrapper( viewKeyStruct::receiverIsLocalString(), &m_receiverIsLocal ). - setInputFlag( InputFlags::FALSE ). - setSizedFromParent( 0 ). - setDescription( "Flag that indicates whether the receiver is local to this MPI rank" ); - registerWrapper( viewKeyStruct::displacementXNp1AtReceiversString(), &m_displacementXNp1AtReceivers ). setInputFlag( InputFlags::FALSE ). setSizedFromParent( 0 ). @@ -268,9 +242,9 @@ void ElasticWaveEquationSEM::postProcessInput() localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 ); m_receiverIsLocal.resize( numReceiversGlobal ); - m_displacementXNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); - m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); + m_displacementXNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); + m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); m_sourceValue.resize( nsamples, numSourcesGlobal ); /// The receivers are initialized to zero. @@ -435,26 +409,9 @@ void ElasticWaveEquationSEM::computeDAS ( arrayView2d< real32 > const xCompRcv, } ); } - /// temporary output to txt - if( m_outputSeismoTrace == 1 ) - { - forAll< serialPolicy >( numReceiversGlobal, [=] ( localIndex const ircv ) - { - if( receiverIsLocal[ircv] == 1 ) - { - std::ofstream f( GEOS_FMT( "dasTraceReceiver{:03}.txt", ircv ), std::ios::app ); - for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample ) - { - f<< iSample << " " << zCompRcv[iSample][ircv] << std::endl; - } - f.close(); - } - } ); - } - /// resize the receiver arrays by dropping the extra pair to avoid confusion /// the remaining x-component contains DAS data, the other components are set to zero - m_displacementXNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); + m_displacementXNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); arrayView2d< real32 > const dasReceiver = m_displacementXNp1AtReceivers.toView(); forAll< EXEC_POLICY >( numReceiversGlobal, [=] GEOS_HOST_DEVICE ( localIndex const ircv ) { @@ -469,9 +426,9 @@ void ElasticWaveEquationSEM::computeDAS ( arrayView2d< real32 > const xCompRcv, } } ); /// set the y and z components to zero to avoid any confusion - m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); + m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); m_displacementYNp1AtReceivers.zero(); - m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); + m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); m_displacementZNp1AtReceivers.zero(); } @@ -587,6 +544,8 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() } ); } ); + WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal ); + WaveSolverUtils::initTrace( "dasTraceReceiver", getName(), m_linearDASGeometry.size( 0 ), m_receiverIsLocal ); } @@ -694,8 +653,6 @@ real64 ElasticWaveEquationSEM::explicitStepInternal( real64 const & time_n, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time_n, dt, cycleNumber ); - GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -780,14 +737,16 @@ real64 ElasticWaveEquationSEM::explicitStepInternal( real64 const & time_n, true ); // compute the seismic traces since last step. - arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView(); - arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView(); - arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView(); + arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView(); + arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView(); + arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView(); computeAllSeismoTraces( time_n, dt, ux_np1, ux_n, uXReceivers ); computeAllSeismoTraces( time_n, dt, uy_np1, uy_n, uYReceivers ); computeAllSeismoTraces( time_n, dt, uz_np1, uz_n, uZReceivers ); + incrementIndexSeismoTrace( time_n ); + forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) { ux_nm1[a] = ux_n[a]; @@ -804,14 +763,8 @@ real64 ElasticWaveEquationSEM::explicitStepInternal( real64 const & time_n, rhsy[a] = 0.0; rhsz[a] = 0.0; } ); - - // increment m_indexSeismoTrace - while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) - { - m_indexSeismoTrace++; - } - } ); + return dt; } @@ -831,50 +784,32 @@ void ElasticWaveEquationSEM::cleanup( real64 const time_n, arrayView1d< string const > const & ) { NodeManager & nodeManager = mesh.getNodeManager(); - arrayView1d< real32 const > const ux_n = nodeManager.getField< fields::Displacementx_n >(); + arrayView1d< real32 const > const ux_n = nodeManager.getField< fields::Displacementx_n >(); arrayView1d< real32 const > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >(); - arrayView1d< real32 const > const uy_n = nodeManager.getField< fields::Displacementy_n >(); + arrayView1d< real32 const > const uy_n = nodeManager.getField< fields::Displacementy_n >(); arrayView1d< real32 const > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >(); - arrayView1d< real32 const > const uz_n = nodeManager.getField< fields::Displacementz_n >(); + arrayView1d< real32 const > const uz_n = nodeManager.getField< fields::Displacementz_n >(); arrayView1d< real32 const > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >(); - arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView(); - arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView(); - arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView(); + arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView(); + arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView(); + arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView(); + + computeAllSeismoTraces( time_n, 0.0, ux_np1, ux_n, uXReceivers ); + computeAllSeismoTraces( time_n, 0.0, uy_np1, uy_n, uYReceivers ); + computeAllSeismoTraces( time_n, 0.0, uz_np1, uz_n, uZReceivers ); - computeAllSeismoTraces( time_n, 0, ux_np1, ux_n, uXReceivers ); - computeAllSeismoTraces( time_n, 0, uy_np1, uy_n, uYReceivers ); - computeAllSeismoTraces( time_n, 0, uz_np1, uz_n, uZReceivers ); + WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, uXReceivers, uYReceivers, uZReceivers ); /// Compute DAS data if requested /// Pairs of receivers are assumed to be modeled ( see WaveSolverBase::initializeDAS() ) if( m_useDAS ) { computeDAS( uXReceivers, uYReceivers, uZReceivers ); + WaveSolverUtils::writeSeismoTrace( "dasTraceReceiver", getName(), m_outputSeismoTrace, m_linearDASGeometry.size( 0 ), + m_receiverIsLocal, m_nsamplesSeismoTrace, uZReceivers ); } } ); - - // increment m_indexSeismoTrace - while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) - { - m_indexSeismoTrace++; - } -} - -void ElasticWaveEquationSEM::computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ) -{ - localIndex indexSeismoTrace = m_indexSeismoTrace; - - for( real64 timeSeismo; - (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace; - indexSeismoTrace++ ) - { - WaveSolverUtils::computeSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal, - m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers ); - } } void ElasticWaveEquationSEM::initializePML() diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp index 15f12013b93..7fd82f1dc56 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp @@ -87,21 +87,6 @@ class ElasticWaveEquationSEM : public WaveSolverBase */ void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhsx, arrayView1d< real32 > const rhsy, arrayView1d< real32 > const rhsz ); - /** - * TODO: move implementation into WaveSolverBase - * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt - * @param time_n the time corresponding to the field values at iteration n - * @param dt the simulation timestep - * @param var_np1 the field values at time_n + dt - * @param var_n the field values at time_n - * @param varAtreceivers the array holding the trace values, where the output is written - */ - virtual void computeAllSeismoTraces( real64 const time_n, - real64 const dt, - arrayView1d< real32 const > const var_np1, - arrayView1d< real32 const > const var_n, - arrayView2d< real32 > varAtReceivers ); - /** * TODO: move implementation into WaveSolverBase once 'm_receiverIsLocal' is also moved * @brief Compute DAS data from the appropriate three-component receiver pairs @@ -123,16 +108,8 @@ class ElasticWaveEquationSEM : public WaveSolverBase real64 const eventProgress, DomainPartition & domain ) override; - struct viewKeyStruct : SolverBase::viewKeyStruct + struct viewKeyStruct : WaveSolverBase::viewKeyStruct { - static constexpr char const * sourceNodeIdsString() { return "sourceNodeIds"; } - static constexpr char const * sourceConstantsString() { return "sourceConstants"; } - static constexpr char const * sourceIsAccessibleString() { return "sourceIsAccessible"; } - - static constexpr char const * receiverNodeIdsString() { return "receiverNodeIds"; } - static constexpr char const * receiverConstantsString() {return "receiverConstants"; } - static constexpr char const * receiverIsLocalString() { return "receiverIsLocal"; } - static constexpr char const * displacementXNp1AtReceiversString() { return "displacementXNp1AtReceivers"; } static constexpr char const * displacementYNp1AtReceiversString() { return "displacementYNp1AtReceivers"; } static constexpr char const * displacementZNp1AtReceiversString() { return "displacementZNp1AtReceivers"; } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp index 83d42014110..78d6cb1006a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp @@ -65,6 +65,11 @@ WaveSolverBase::WaveSolverBase( const std::string & name, setApplyDefaultValue( -1 ). setDescription( "Source time delay (1 / f0 by default)" ); + registerWrapper( viewKeyStruct::timeSourceFrequencyString(), &m_timeSourceFrequency ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "Central frequency for the time source" ); + registerWrapper( viewKeyStruct::rickerOrderString(), &m_rickerOrder ). setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 2 ). @@ -120,7 +125,6 @@ WaveSolverBase::WaveSolverBase( const std::string & name, setApplyDefaultValue( -80 ). setDescription( "Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory)" ); - registerWrapper( viewKeyStruct::usePMLString(), &m_usePML ). setInputFlag( InputFlags::FALSE ). setApplyDefaultValue( 0 ). @@ -156,7 +160,7 @@ WaveSolverBase::WaveSolverBase( const std::string & name, setSizedFromParent( 0 ). setDescription( "Indices of the nodes (in the right order) for each receiver point" ); - registerWrapper( viewKeyStruct::sourceConstantsString(), &m_sourceConstants ). + registerWrapper( viewKeyStruct::receiverConstantsString(), &m_receiverConstants ). setInputFlag( InputFlags::FALSE ). setSizedFromParent( 0 ). setDescription( "Constant part of the receiver for the nodes listed in m_receiverNodeIds" ); @@ -166,7 +170,15 @@ WaveSolverBase::WaveSolverBase( const std::string & name, setSizedFromParent( 0 ). setDescription( "Flag that indicates whether the receiver is local to this MPI rank" ); + registerWrapper( viewKeyStruct::receiverRegionString(), &m_receiverRegion ). + setInputFlag( InputFlags::FALSE ). + setSizedFromParent( 0 ). + setDescription( "Region containing the receivers" ); + registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ). + setInputFlag( InputFlags::FALSE ). + setSizedFromParent( 0 ). + setDescription( "Element containing the receivers" ); } WaveSolverBase::~WaveSolverBase() @@ -241,38 +253,32 @@ void WaveSolverBase::postProcessInput() m_usePML = counter; if( m_linearDASGeometry.size( 1 ) > 0 ) - { m_useDAS = 1; - } if( m_useDAS ) { GEOS_LOG_LEVEL_RANK_0( 1, "Modeling linear DAS data is activated" ); GEOS_ERROR_IF( m_linearDASGeometry.size( 1 ) != 3, - getWrapperDataContext( viewKeyStruct::linearDASGeometryString() ) << - ": Invalid number of geometry parameters for the linear DAS fiber. Three parameters are required: dip, azimuth, gauge length" ); + "Invalid number of geometry parameters for the linear DAS fiber. Three parameters are required: dip, azimuth, gauge length" ); GEOS_ERROR_IF( m_linearDASGeometry.size( 0 ) != m_receiverCoordinates.size( 0 ), - getWrapperDataContext( viewKeyStruct::linearDASGeometryString() ) << - ": Invalid number of geometry parameters instances for the linear DAS fiber. It should match the number of receivers." ); + "Invalid number of geometry parameters instances for the linear DAS fiber. It should match the number of receivers." ); /// initialize DAS geometry initializeDAS(); } - GEOS_THROW_IF( m_sourceCoordinates.size( 1 ) != 3, - getWrapperDataContext( viewKeyStruct::sourceCoordinatesString() ) << - ": Invalid number of physical coordinates for the sources", + GEOS_THROW_IF( m_sourceCoordinates.size( 0 ) > 0 && m_sourceCoordinates.size( 1 ) != 3, + "Invalid number of physical coordinates for the sources", InputError ); - GEOS_THROW_IF( m_receiverCoordinates.size( 1 ) != 3, - getWrapperDataContext( viewKeyStruct::receiverCoordinatesString() ) << - ": Invalid number of physical coordinates for the receivers", + GEOS_THROW_IF( m_receiverCoordinates.size( 0 ) > 0 && m_receiverCoordinates.size( 1 ) != 3, + "Invalid number of physical coordinates for the receivers", InputError ); - EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" ); + EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() ); real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() ); real64 dt = 0; @@ -285,7 +291,7 @@ void WaveSolverBase::postProcessInput() } } - GEOS_THROW_IF( dt < epsilonLoc*maxTime, getDataContext() << ": Value for dt: " << dt <<" is smaller than local threshold: " << epsilonLoc, std::runtime_error ); + GEOS_THROW_IF( dt < epsilonLoc * maxTime, getDataContext() << ": Value for dt: " << dt <<" is smaller than local threshold: " << epsilonLoc, std::runtime_error ); if( m_dtSeismoTrace > 0 ) { @@ -295,7 +301,7 @@ void WaveSolverBase::postProcessInput() { m_nsamplesSeismoTrace = 0; } - localIndex const nsamples = int( (maxTime-minTime) /dt) + 1; + localIndex const nsamples = int( (maxTime - minTime) / dt) + 1; localIndex const numSourcesGlobal = m_sourceCoordinates.size( 0 ); m_sourceValue.resize( nsamples, numSourcesGlobal ); @@ -396,7 +402,67 @@ localIndex WaveSolverBase::getNumNodesPerElem() } ); return numNodesPerElem; +} +void WaveSolverBase::incrementIndexSeismoTrace( real64 const time_n ) +{ + while( (m_dtSeismoTrace * m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) + { + m_indexSeismoTrace++; + } +} + +void WaveSolverBase::computeAllSeismoTraces( real64 const time_n, + real64 const dt, + arrayView1d< real32 const > const var_np1, + arrayView1d< real32 const > const var_n, + arrayView2d< real32 > varAtReceivers ) +{ + /* + * In forward case we compute seismo if time_n + dt is the first time + * step after the timeSeismo to write. + * + * time_n timeSeismo time_n + dt + * ---|--------------|-------------| + * + * In backward (time_n goes decreasing) case we compute seismo if + * time_n is the last time step before the timeSeismo to write. + * + * time_n - dt timeSeismo time_n + * ---|--------------|-------------| + */ + + if( m_nsamplesSeismoTrace == 0 ) + return; + integer const dir = m_forward ? +1 : -1; + for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ ) + { + real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo); + if( dir * timeSeismo > dir * (time_n + epsilonLoc)) + break; + WaveSolverUtils::computeSeismoTrace( time_n, dir * dt, timeSeismo, iSeismo, m_receiverNodeIds, + m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers ); + } +} + +void WaveSolverBase::compute2dVariableAllSeismoTraces( localIndex const regionIndex, + real64 const time_n, + real64 const dt, + arrayView2d< real32 const > const var_np1, + arrayView2d< real32 const > const var_n, + arrayView2d< real32 > varAtReceivers ) +{ + if( m_nsamplesSeismoTrace == 0 ) + return; + integer const dir = m_forward ? +1 : -1; + for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ ) + { + real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo); + if( dir * timeSeismo > dir * (time_n + epsilonLoc)) + break; + WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dir * dt, regionIndex, m_receiverRegion, timeSeismo, iSeismo, m_rcvElem, + m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers ); + } } bool WaveSolverBase::directoryExists( std::string const & directoryName ) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp index f4831b8499c..820f5e8f58c 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp @@ -27,7 +27,7 @@ #if !defined( GEOS_USE_HIP ) #include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" #endif - +#include "WaveSolverUtils.hpp" #if !defined( GEOS_USE_HIP ) #define SEM_FE_TYPES \ @@ -49,8 +49,9 @@ class WaveSolverBase : public SolverBase { public: - using EXEC_POLICY = parallelDevicePolicy< >; - using wsCoordType = real32; + static constexpr real64 epsilonLoc = WaveSolverUtils::epsilonLoc; + using EXEC_POLICY = WaveSolverUtils::EXEC_POLICY; + using wsCoordType = WaveSolverUtils::wsCoordType; WaveSolverBase( const std::string & name, Group * const parent ); @@ -113,14 +114,11 @@ class WaveSolverBase : public SolverBase static constexpr char const * usePMLString() { return "usePML"; } static constexpr char const * parametersPMLString() { return "parametersPML"; } + static constexpr char const * receiverElemString() { return "rcvElem"; } + static constexpr char const * receiverRegionString() { return "receiverRegion"; } static constexpr char const * freeSurfaceString() { return "FreeSurface"; } }; - /** - * @brief Safeguard for timeStep. Used to avoid memory issue due to too small value. - */ - static constexpr real64 epsilonLoc = 1e-8; - /** * @brief Re-initialize source and receivers positions in the mesh, and resize the pressureNp1_at_receivers array */ @@ -154,6 +152,36 @@ class WaveSolverBase : public SolverBase */ virtual void initializePML() = 0; + virtual void incrementIndexSeismoTrace( real64 const time_n ); + + /** + * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt + * @param time_n the time corresponding to the field values pressure_n + * @param dt the simulation timestep + * @param var_np1 the field values at time_n + dt + * @param var_n the field values at time_n + * @param varAtreceivers the array holding the trace values, where the output is written + */ + virtual void computeAllSeismoTraces( real64 const time_n, + real64 const dt, + arrayView1d< real32 const > const var_np1, + arrayView1d< real32 const > const var_n, + arrayView2d< real32 > varAtReceivers ); + /** + * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt + * @param time_n the time corresponding to the field values pressure_n + * @param dt the simulation timestep + * @param var_np1 the field values at time_n + dt + * @param var_n the field values at time_n + * @param varAtreceivers the array holding the trace values, where the output is written + */ + virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex, + real64 const time_n, + real64 const dt, + arrayView2d< real32 const > const var_np1, + arrayView2d< real32 const > const var_n, + arrayView2d< real32 > varAtReceivers ); + /** * @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 @@ -161,8 +189,6 @@ class WaveSolverBase : public SolverBase */ virtual void applyPML( real64 const time, DomainPartition & domain ) = 0; - - /** * @brief Locate sources and receivers positions in the mesh elements, evaluate the basis functions at each point and save them to the * corresponding elements nodes. @@ -202,7 +228,6 @@ class WaveSolverBase : public SolverBase virtual void registerDataOnMesh( Group & meshBodies ) override; - localIndex getNumNodesPerElem(); /// Coordinates of the sources in the mesh @@ -224,7 +249,7 @@ class WaveSolverBase : public SolverBase localIndex m_rickerOrder; /// Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise - localIndex m_outputSeismoTrace; + integer m_outputSeismoTrace; /// Time step for seismoTrace output real64 m_dtSeismoTrace; @@ -271,6 +296,12 @@ class WaveSolverBase : public SolverBase /// Flag that indicates whether the receiver is local or not to the MPI rank array1d< localIndex > m_receiverIsLocal; + /// Array containing the elements which contain a receiver + array1d< localIndex > m_rcvElem; + + /// Array containing the elements which contain the region which the receiver belongs + array1d< localIndex > m_receiverRegion; + /// Flag to enable LIFO localIndex m_enableLifo; @@ -310,7 +341,7 @@ class WaveSolverBase : public SolverBase namespace fields { -using reference32Type = array2d< WaveSolverBase::wsCoordType, nodes::REFERENCE_POSITION_PERM >; +using reference32Type = array2d< WaveSolverUtils::wsCoordType, nodes::REFERENCE_POSITION_PERM >; DECLARE_FIELD( referencePosition32, "referencePosition32", reference32Type, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp index 2915dfb24fa..ca52e3af4c0 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp @@ -20,14 +20,21 @@ #ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERUTILS_HPP_ #define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERUTILS_HPP_ +#include "mesh/utilities/ComputationalGeometry.hpp" #include "fileIO/Outputs/OutputBase.hpp" -#include "WaveSolverBase.hpp" +#include "LvArray/src/tensorOps.hpp" namespace geos { struct WaveSolverUtils { + static constexpr real64 epsilonLoc = 1e-8; + static constexpr real64 eps64 = std::numeric_limits< real64 >::epsilon(); + static constexpr real32 eps32 = std::numeric_limits< real32 >::epsilon(); + + using EXEC_POLICY = parallelDevicePolicy< >; + using wsCoordType = real32; GEOS_HOST_DEVICE static real32 evaluateRicker( real64 const time_n, real32 const f0, real32 const t0, localIndex const order ) @@ -68,78 +75,124 @@ struct WaveSolverUtils return pulse; } - static void writeSeismoTrace( localIndex iSeismo, - arrayView2d< real64 const > const receiverConstants, + /** + * @brief Initialize (clear) the trace file. + */ + static void initTrace( char const * prefix, + string const & name, + localIndex const nReceivers, + arrayView1d< localIndex const > const receiverIsLocal ) + { + string const outputDir = OutputBase::getOutputDirectory(); + RAJA::ReduceSum< ReducePolicy< serialPolicy >, localIndex > count( 0 ); + + forAll< serialPolicy >( nReceivers, [=] ( localIndex const ircv ) + { + if( receiverIsLocal[ircv] == 1 ) + { + count += 1; + string const fn = joinPath( outputDir, GEOS_FMT( "{}_{}_{:03}.txt", prefix, name, ircv ) ); + std::ofstream f( fn, std::ios::out | std::ios::trunc ); + } + } ); + + localIndex const total = MpiWrapper::sum( count.get() ); + GEOS_ERROR_IF( nReceivers != total, GEOS_FMT( ": Invalid distribution of receivers: nReceivers={} != MPI::sum={}.", nReceivers, total ) ); + } + + /** + * @brief Convenient helper for 3D vectors calling 3 times the scalar version. + */ + static void writeSeismoTraceVector( char const * prefix, + string const & name, + bool const outputSeismoTrace, + localIndex const nReceivers, + arrayView1d< localIndex const > const receiverIsLocal, + localIndex const nsamplesSeismoTrace, + arrayView2d< real32 const > const varAtReceiversx, + arrayView2d< real32 const > const varAtReceiversy, + arrayView2d< real32 const > const varAtReceiversz ) + { + writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversx ); + writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversy ); + writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversz ); + } + + /** + * @brief Write the seismo traces to a file. + */ + static void writeSeismoTrace( char const * prefix, + string const & name, + bool const outputSeismoTrace, + localIndex const nReceivers, arrayView1d< localIndex const > const receiverIsLocal, localIndex const nsamplesSeismoTrace, - localIndex const outputSeismoTrace, - arrayView2d< real32 > varAtReceivers ) + arrayView2d< real32 const > const varAtReceivers ) { - if( iSeismo == nsamplesSeismoTrace - 1 ) + if( !outputSeismoTrace ) return; + + string const outputDir = OutputBase::getOutputDirectory(); + forAll< serialPolicy >( nReceivers, [=] ( localIndex const ircv ) { - string const outputDir = OutputBase::getOutputDirectory(); - forAll< serialPolicy >( receiverConstants.size( 0 ), [=] ( localIndex const ircv ) + if( receiverIsLocal[ircv] == 1 ) { - if( outputSeismoTrace == 1 ) + string const fn = joinPath( outputDir, GEOS_FMT( "{}_{}_{:03}.txt", prefix, name, ircv ) ); + std::ofstream f( fn, std::ios::app ); + if( f ) { - if( receiverIsLocal[ircv] == 1 ) + for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample ) { - string const fn = joinPath( outputDir, GEOS_FMT( "seismoTraceReceiver{:03}.txt", ircv ) ); - std::ofstream f( fn, std::ios::app ); - if( !f ) - { - GEOS_WARNING( GEOS_FMT( "Failed to open output file {}", fn ) ); - return; - } - for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample ) - { - f << iSample << " " << varAtReceivers[iSample][ircv] << std::endl; - } - f.close(); + // index - time - value + f << iSample << " " << varAtReceivers[iSample][nReceivers] << " " << varAtReceivers[iSample][ircv] << std::endl; } + f.close(); } - } ); - } + else + { + GEOS_WARNING( GEOS_FMT( "Failed to open output file {}", fn ) ); + } + } + } ); } + /** + * @brief Compute the seismo traces. + */ static void computeSeismoTrace( real64 const time_n, real64 const dt, real64 const timeSeismo, - localIndex iSeismo, + localIndex const iSeismo, arrayView2d< localIndex const > const receiverNodeIds, arrayView2d< real64 const > const receiverConstants, arrayView1d< localIndex const > const receiverIsLocal, - localIndex const nsamplesSeismoTrace, - localIndex const outputSeismoTrace, arrayView1d< real32 const > const var_np1, arrayView1d< real32 const > const var_n, arrayView2d< real32 > varAtReceivers ) { real64 const time_np1 = time_n + dt; - real32 const a1 = (LvArray::math::abs( dt ) < WaveSolverBase::epsilonLoc ) ? 1.0 : (time_np1 - timeSeismo)/dt; + real32 const a1 = abs( dt ) < epsilonLoc ? 1.0 : (time_np1 - timeSeismo) / dt; real32 const a2 = 1.0 - a1; - if( nsamplesSeismoTrace > 0 ) + localIndex const nReceivers = receiverConstants.size( 0 ); + + forAll< EXEC_POLICY >( nReceivers, [=] GEOS_HOST_DEVICE ( localIndex const ircv ) { - forAll< WaveSolverBase::EXEC_POLICY >( receiverConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const ircv ) + if( receiverIsLocal[ircv] == 1 ) { - if( receiverIsLocal[ircv] == 1 ) + real32 vtmp_np1 = 0.0, vtmp_n = 0.0; + for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode ) { - varAtReceivers[iSeismo][ircv] = 0.0; - real32 vtmp_np1 = 0.0, vtmp_n = 0.0; - for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode ) - { - vtmp_np1 += var_np1[receiverNodeIds[ircv][inode]] * receiverConstants[ircv][inode]; - vtmp_n += var_n[receiverNodeIds[ircv][inode]] * receiverConstants[ircv][inode]; - } - // linear interpolation between the pressure value at time_n and time_(n+1) - varAtReceivers[iSeismo][ircv] = a1*vtmp_n + a2*vtmp_np1; + vtmp_np1 += var_np1[receiverNodeIds( ircv, inode )] * receiverConstants( ircv, inode ); + vtmp_n += var_n[receiverNodeIds( ircv, inode )] * receiverConstants( ircv, inode ); } - } ); - } - - writeSeismoTrace( iSeismo, receiverConstants, receiverIsLocal, nsamplesSeismoTrace, outputSeismoTrace, varAtReceivers ); + // linear interpolation between the pressure value at time_n and time_{n+1} + varAtReceivers( iSeismo, ircv ) = a1 * vtmp_n + a2 * vtmp_np1; + // NOTE: varAtReceivers has size(1) = numReceiversGlobal + 1, this does not OOB + // left in the forAll loop for sync issues since the following does not depend on `ircv` + varAtReceivers( iSeismo, nReceivers ) = a1 * time_n + a2 * time_np1; + } + } ); } static void compute2dVariableSeismoTrace( real64 const time_n, @@ -147,50 +200,43 @@ struct WaveSolverUtils localIndex const regionIndex, arrayView1d< localIndex const > const receiverRegion, real64 const timeSeismo, - localIndex iSeismo, + localIndex const iSeismo, arrayView1d< localIndex const > const rcvElem, arrayView2d< real64 const > const receiverConstants, arrayView1d< localIndex const > const receiverIsLocal, - localIndex const nsamplesSeismoTrace, - localIndex const outputSeismoTrace, arrayView2d< real32 const > const var_np1, arrayView2d< real32 const > const var_n, arrayView2d< real32 > varAtReceivers ) { real64 const time_np1 = time_n+dt; - real32 const a1 = (dt < WaveSolverBase::epsilonLoc) ? 1.0 : (time_np1 - timeSeismo)/dt; + real32 const a1 = dt < epsilonLoc ? 1.0 : (time_np1 - timeSeismo) / dt; real32 const a2 = 1.0 - a1; - if( nsamplesSeismoTrace > 0 ) + localIndex const nReceivers = receiverConstants.size( 0 ); + + forAll< EXEC_POLICY >( nReceivers, [=] GEOS_HOST_DEVICE ( localIndex const ircv ) { - forAll< WaveSolverBase::EXEC_POLICY >( receiverConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const ircv ) + if( receiverIsLocal[ircv] == 1 ) { - if( receiverIsLocal[ircv] == 1 ) + if( receiverRegion[ircv] == regionIndex ) { - if( receiverRegion[ircv] == regionIndex ) + real32 vtmp_np1 = 0.0, vtmp_n = 0.0; + for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode ) { - varAtReceivers[iSeismo][ircv] = 0.0; - real32 vtmp_np1 = 0.0, vtmp_n = 0.0; - for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode ) - { - vtmp_np1 += var_np1[rcvElem[ircv]][inode] * receiverConstants[ircv][inode]; - vtmp_n += var_n[rcvElem[ircv]][inode] * receiverConstants[ircv][inode]; - } - // linear interpolation between the pressure value at time_n and time_(n+1) - varAtReceivers[iSeismo][ircv] = a1*vtmp_n + a2*vtmp_np1; + vtmp_np1 += var_np1( rcvElem[ircv], inode ) * receiverConstants( ircv, inode ); + vtmp_n += var_n( rcvElem[ircv], inode ) * receiverConstants( ircv, inode ); } + // linear interpolation between the pressure value at time_n and time_{n+1} + varAtReceivers( iSeismo, ircv ) = a1 * vtmp_n + a2 * vtmp_np1; + // NOTE: varAtReceivers has size(1) = numReceiversGlobal + 1, this does not OOB + // left in the forAll loop for sync issues since the following does not depend on `ircv` + varAtReceivers( iSeismo, nReceivers ) = a1 * time_n + a2 * time_np1; } - } ); - } - - writeSeismoTrace( iSeismo, receiverConstants, receiverIsLocal, nsamplesSeismoTrace, outputSeismoTrace, varAtReceivers ); + } + } ); } - /** - * @brief Check if the source point is inside an element or not - */ - /** * @brief Check if the source point is inside an element or not * @param numFacesPerElem number of face on an element @@ -201,7 +247,6 @@ struct WaveSolverUtils * @param coords coordinate of the point * @return true if coords is inside the element */ - GEOS_HOST_DEVICE static bool locateSourceElement( real64 const numFacesPerElem, @@ -237,35 +282,32 @@ struct WaveSolverUtils localIndex const s = computationalGeometry::sign( LvArray::tensorOps::AiBi< 3 >( faceNormalOnFace, faceCenterOnFace )); // all dot products should be non-negative (we enforce outward normals) - if( s < 0 ) - { - return false; - } + if( s < 0 ) return false; } return true; } -/** - * @brief Convert a mesh element point coordinate into a coordinate on the reference element - * @tparam FE_TYPE finite element type - * @param[in] coords coordinate of the point - * @param[in] elemsToNodes map to obtaint global nodes from element index - * @param[in] X array of mesh nodes coordinates - * @param[out] coordsOnRefElem to contain the coordinate computed in the reference element - */ + /** + * @brief Convert a mesh element point coordinate into a coordinate on the reference element + * @tparam FE_TYPE finite element type + * @param[in] coords coordinate of the point + * @param[in] elemsToNodes map to obtaint global nodes from element index + * @param[in] nodeCoords array of mesh nodes coordinates + * @param[out] coordsOnRefElem to contain the coordinate computed in the reference element + */ template< typename FE_TYPE > GEOS_HOST_DEVICE static void computeCoordinatesOnReferenceElement( real64 const (&coords)[3], arraySlice1d< localIndex const, cells::NODE_MAP_USD - 1 > const elemsToNodes, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, real64 (& coordsOnRefElem)[3] ) { real64 xLocal[FE_TYPE::numNodes][3]{}; for( localIndex a = 0; a < FE_TYPE::numNodes; ++a ) { - LvArray::tensorOps::copy< 3 >( xLocal[a], X[ elemsToNodes[a] ] ); + LvArray::tensorOps::copy< 3 >( xLocal[a], nodeCoords[ elemsToNodes[a] ] ); } // coordsOnRefElem = invJ*(coords-coordsNode_0) real64 invJ[3][3]{}; diff --git a/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst b/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst index 9ad3b2e0db1..a5100a7e2f7 100644 --- a/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst +++ b/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst @@ -23,7 +23,7 @@ shotIndex integer 0 Set the current shot for tem sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. timeSourceDelay real32 -1 Source time delay (1 / f0 by default) -timeSourceFrequency real32 required Central frequency for the time source +timeSourceFrequency real32 0 Central frequency for the time source LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` ========================= ============== ========== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/AcousticFirstOrderSEM_other.rst b/src/coreComponents/schema/docs/AcousticFirstOrderSEM_other.rst index 37156e3e38a..24bd118a9d9 100644 --- a/src/coreComponents/schema/docs/AcousticFirstOrderSEM_other.rst +++ b/src/coreComponents/schema/docs/AcousticFirstOrderSEM_other.rst @@ -8,10 +8,11 @@ maxStableDt real64 meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. pressureNp1AtReceivers real32_array2d Pressure value at each receiver for each timestep rcvElem integer_array Element containing the receivers +receiverConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point receiverRegion integer_array Region containing the receivers -sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds +sourceConstants real64_array2d Constant part of the source for the nodes listed in m_sourceNodeIds sourceElem integer_array Element containing the sources sourceIsAccessible integer_array Flag that indicates whether the source is local to this MPI rank sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point diff --git a/src/coreComponents/schema/docs/AcousticSEM.rst b/src/coreComponents/schema/docs/AcousticSEM.rst index 9ad3b2e0db1..a5100a7e2f7 100644 --- a/src/coreComponents/schema/docs/AcousticSEM.rst +++ b/src/coreComponents/schema/docs/AcousticSEM.rst @@ -23,7 +23,7 @@ shotIndex integer 0 Set the current shot for tem sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. timeSourceDelay real32 -1 Source time delay (1 / f0 by default) -timeSourceFrequency real32 required Central frequency for the time source +timeSourceFrequency real32 0 Central frequency for the time source LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` ========================= ============== ========== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/AcousticSEM_other.rst b/src/coreComponents/schema/docs/AcousticSEM_other.rst index a5f09a365c9..0f439ecdf9b 100644 --- a/src/coreComponents/schema/docs/AcousticSEM_other.rst +++ b/src/coreComponents/schema/docs/AcousticSEM_other.rst @@ -7,9 +7,12 @@ indexSeismoTrace integer maxStableDt real64 Value of the Maximum Stable Timestep for this solver. meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. pressureNp1AtReceivers real32_array2d Pressure value at each receiver for each timestep +rcvElem integer_array Element containing the receivers +receiverConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point -sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds +receiverRegion integer_array Region containing the receivers +sourceConstants real64_array2d Constant part of the source for the nodes listed in m_sourceNodeIds sourceIsAccessible integer_array Flag that indicates whether the source is local to this MPI rank sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point sourceValue real32_array2d Source Value of the sources diff --git a/src/coreComponents/schema/docs/AcousticVTISEM.rst b/src/coreComponents/schema/docs/AcousticVTISEM.rst index 9ad3b2e0db1..a5100a7e2f7 100644 --- a/src/coreComponents/schema/docs/AcousticVTISEM.rst +++ b/src/coreComponents/schema/docs/AcousticVTISEM.rst @@ -23,7 +23,7 @@ shotIndex integer 0 Set the current shot for tem sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. timeSourceDelay real32 -1 Source time delay (1 / f0 by default) -timeSourceFrequency real32 required Central frequency for the time source +timeSourceFrequency real32 0 Central frequency for the time source LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` ========================= ============== ========== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/AcousticVTISEM_other.rst b/src/coreComponents/schema/docs/AcousticVTISEM_other.rst index a5f09a365c9..0f439ecdf9b 100644 --- a/src/coreComponents/schema/docs/AcousticVTISEM_other.rst +++ b/src/coreComponents/schema/docs/AcousticVTISEM_other.rst @@ -7,9 +7,12 @@ indexSeismoTrace integer maxStableDt real64 Value of the Maximum Stable Timestep for this solver. meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. pressureNp1AtReceivers real32_array2d Pressure value at each receiver for each timestep +rcvElem integer_array Element containing the receivers +receiverConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point -sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds +receiverRegion integer_array Region containing the receivers +sourceConstants real64_array2d Constant part of the source for the nodes listed in m_sourceNodeIds sourceIsAccessible integer_array Flag that indicates whether the source is local to this MPI rank sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point sourceValue real32_array2d Source Value of the sources diff --git a/src/coreComponents/schema/docs/BlackOilFluid.rst b/src/coreComponents/schema/docs/BlackOilFluid.rst index c84a4c0613d..75b10845be1 100644 --- a/src/coreComponents/schema/docs/BlackOilFluid.rst +++ b/src/coreComponents/schema/docs/BlackOilFluid.rst @@ -3,6 +3,7 @@ ======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== Name Type Default Description ======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. componentMolarWeight real64_array required Component molar weights componentNames string_array {} List of component names hydrocarbonFormationVolFactorTableNames string_array {} | List of formation volume factor TableFunction names from the Functions block. diff --git a/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid.rst b/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid.rst index 2010eb6efc8..c4a713166eb 100644 --- a/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid.rst +++ b/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid.rst @@ -1,14 +1,15 @@ -==================== ============ ======== ============================================================================== -Name Type Default Description -==================== ============ ======== ============================================================================== -componentMolarWeight real64_array {0} Component molar weights -componentNames string_array {} List of component names -flashModelParaFile path required Name of the file defining the parameters of the flash model -name string required A name is required for any non-unique nodes -phaseNames string_array {} List of fluid phases -phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models -==================== ============ ======== ============================================================================== +==================== ============ ======== ============================================================================================================ +Name Type Default Description +==================== ============ ======== ============================================================================================================ +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. +componentMolarWeight real64_array {0} Component molar weights +componentNames string_array {} List of component names +flashModelParaFile path required Name of the file defining the parameters of the flash model +name string required A name is required for any non-unique nodes +phaseNames string_array {} List of fluid phases +phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models +==================== ============ ======== ============================================================================================================ diff --git a/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid.rst b/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid.rst index 2010eb6efc8..c4a713166eb 100644 --- a/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid.rst +++ b/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid.rst @@ -1,14 +1,15 @@ -==================== ============ ======== ============================================================================== -Name Type Default Description -==================== ============ ======== ============================================================================== -componentMolarWeight real64_array {0} Component molar weights -componentNames string_array {} List of component names -flashModelParaFile path required Name of the file defining the parameters of the flash model -name string required A name is required for any non-unique nodes -phaseNames string_array {} List of fluid phases -phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models -==================== ============ ======== ============================================================================== +==================== ============ ======== ============================================================================================================ +Name Type Default Description +==================== ============ ======== ============================================================================================================ +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. +componentMolarWeight real64_array {0} Component molar weights +componentNames string_array {} List of component names +flashModelParaFile path required Name of the file defining the parameters of the flash model +name string required A name is required for any non-unique nodes +phaseNames string_array {} List of fluid phases +phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models +==================== ============ ======== ============================================================================================================ diff --git a/src/coreComponents/schema/docs/CO2BrinePhillipsFluid.rst b/src/coreComponents/schema/docs/CO2BrinePhillipsFluid.rst index 2010eb6efc8..c4a713166eb 100644 --- a/src/coreComponents/schema/docs/CO2BrinePhillipsFluid.rst +++ b/src/coreComponents/schema/docs/CO2BrinePhillipsFluid.rst @@ -1,14 +1,15 @@ -==================== ============ ======== ============================================================================== -Name Type Default Description -==================== ============ ======== ============================================================================== -componentMolarWeight real64_array {0} Component molar weights -componentNames string_array {} List of component names -flashModelParaFile path required Name of the file defining the parameters of the flash model -name string required A name is required for any non-unique nodes -phaseNames string_array {} List of fluid phases -phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models -==================== ============ ======== ============================================================================== +==================== ============ ======== ============================================================================================================ +Name Type Default Description +==================== ============ ======== ============================================================================================================ +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. +componentMolarWeight real64_array {0} Component molar weights +componentNames string_array {} List of component names +flashModelParaFile path required Name of the file defining the parameters of the flash model +name string required A name is required for any non-unique nodes +phaseNames string_array {} List of fluid phases +phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models +==================== ============ ======== ============================================================================================================ diff --git a/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid.rst b/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid.rst index 2010eb6efc8..c4a713166eb 100644 --- a/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid.rst +++ b/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid.rst @@ -1,14 +1,15 @@ -==================== ============ ======== ============================================================================== -Name Type Default Description -==================== ============ ======== ============================================================================== -componentMolarWeight real64_array {0} Component molar weights -componentNames string_array {} List of component names -flashModelParaFile path required Name of the file defining the parameters of the flash model -name string required A name is required for any non-unique nodes -phaseNames string_array {} List of fluid phases -phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models -==================== ============ ======== ============================================================================== +==================== ============ ======== ============================================================================================================ +Name Type Default Description +==================== ============ ======== ============================================================================================================ +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. +componentMolarWeight real64_array {0} Component molar weights +componentNames string_array {} List of component names +flashModelParaFile path required Name of the file defining the parameters of the flash model +name string required A name is required for any non-unique nodes +phaseNames string_array {} List of fluid phases +phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models +==================== ============ ======== ============================================================================================================ diff --git a/src/coreComponents/schema/docs/CellElementRegion.rst b/src/coreComponents/schema/docs/CellElementRegion.rst index eefbb5a2c43..ae17d10db84 100644 --- a/src/coreComponents/schema/docs/CellElementRegion.rst +++ b/src/coreComponents/schema/docs/CellElementRegion.rst @@ -3,7 +3,7 @@ =============== ============ ======== =========================================== Name Type Default Description =============== ============ ======== =========================================== -cellBlocks string_array {} (no description available) +cellBlocks string_array required (no description available) coarseningRatio real64 0 (no description available) materialList string_array required List of materials present in this region meshBody string Mesh body that contains this region diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst index c89bf1c12bb..25e64a02616 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst @@ -12,6 +12,7 @@ logLevel integer maxCompFractionChange real64 0.5 Maximum (absolute) change in a component fraction in a Newton iteration maxRelativePressureChange real64 0.5 Maximum (relative) change in pressure in a Newton iteration (expected value between 0 and 1) maxRelativeTemperatureChange real64 0.5 Maximum (relative) change in temperature in a Newton iteration (expected value between 0 and 1) +minCompDens real64 1e-10 Minimum allowed global component density name string required A name is required for any non-unique nodes scalingType geos_CompositionalMultiphaseFVM_ScalingType Global | Solution scaling type.Valid options: | * Global @@ -23,6 +24,8 @@ targetRelativePressureChangeInTimeStep real64 targetRelativeTemperatureChangeInTimeStep real64 0.2 Target (relative) change in temperature in a time step (expected value between 0 and 1) temperature real64 required Temperature useMass integer 0 Use mass formulation instead of molar +useSimpleAccumulation integer 0 Flag indicating whether simple accumulation form is used +useTotalMassEquation integer 1 Flag indicating whether total mass equation is used LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` ========================================= =========================================== ======== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseFluid.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseFluid.rst index cdfebc46753..97e11ca2bb5 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseFluid.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseFluid.rst @@ -1,18 +1,19 @@ -============================ ============== ======== ============================================== -Name Type Default Description -============================ ============== ======== ============================================== -componentAcentricFactor real64_array required Component acentric factors -componentBinaryCoeff real64_array2d {{0}} Table of binary interaction coefficients -componentCriticalPressure real64_array required Component critical pressures -componentCriticalTemperature real64_array required Component critical temperatures -componentMolarWeight real64_array required Component molar weights -componentNames string_array required List of component names -componentVolumeShift real64_array {0} Component volume shifts -equationsOfState string_array required List of equation of state types for each phase -name string required A name is required for any non-unique nodes -phaseNames string_array required List of fluid phases -============================ ============== ======== ============================================== +============================ ============== ======== ============================================================================================================ +Name Type Default Description +============================ ============== ======== ============================================================================================================ +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. +componentAcentricFactor real64_array required Component acentric factors +componentBinaryCoeff real64_array2d {{0}} Table of binary interaction coefficients +componentCriticalPressure real64_array required Component critical pressures +componentCriticalTemperature real64_array required Component critical temperatures +componentMolarWeight real64_array required Component molar weights +componentNames string_array required List of component names +componentVolumeShift real64_array {0} Component volume shifts +equationsOfState string_array required List of equation of state types for each phase +name string required A name is required for any non-unique nodes +phaseNames string_array required List of fluid phases +============================ ============== ======== ============================================================================================================ diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst index 6648f0b0d41..dba62a0e6ab 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst @@ -12,6 +12,7 @@ logLevel integer 0 Log level maxCompFractionChange real64 0.5 Maximum (absolute) change in a component fraction in a Newton iteration maxRelativePressureChange real64 0.5 Maximum (relative) change in pressure in a Newton iteration (expected value between 0 and 1) maxRelativeTemperatureChange real64 0.5 Maximum (relative) change in temperature in a Newton iteration (expected value between 0 and 1) +minCompDens real64 1e-10 Minimum allowed global component density name string required A name is required for any non-unique nodes solutionChangeScalingFactor real64 0.5 Damping factor for solution change targets targetPhaseVolFractionChangeInTimeStep real64 0.2 Target (absolute) change in phase volume fraction in a time step @@ -20,6 +21,8 @@ targetRelativePressureChangeInTimeStep real64 0.2 Target (relative targetRelativeTemperatureChangeInTimeStep real64 0.2 Target (relative) change in temperature in a time step (expected value between 0 and 1) temperature real64 required Temperature useMass integer 0 Use mass formulation instead of molar +useSimpleAccumulation integer 0 Flag indicating whether simple accumulation form is used +useTotalMassEquation integer 1 Flag indicating whether total mass equation is used LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` ========================================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseWell.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseWell.rst index 6af75cdcb4e..a3540906c9a 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseWell.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseWell.rst @@ -11,7 +11,7 @@ maxCompFractionChange real64 1 Maximum (absolute) change in maxRelativePressureChange real64 1 Maximum (relative) change in pressure between two Newton iterations (recommended with rate control) name string required A name is required for any non-unique nodes targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. -useMass integer 0 Use mass formulation instead of molar +useMass integer 0 Use total mass equation LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` WellControls node :ref:`XML_WellControls` diff --git a/src/coreComponents/schema/docs/DeadOilFluid.rst b/src/coreComponents/schema/docs/DeadOilFluid.rst index c84a4c0613d..75b10845be1 100644 --- a/src/coreComponents/schema/docs/DeadOilFluid.rst +++ b/src/coreComponents/schema/docs/DeadOilFluid.rst @@ -3,6 +3,7 @@ ======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== Name Type Default Description ======================================= ============ ======== ===================================================================================================================================================================================================================================================================================================== +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. componentMolarWeight real64_array required Component molar weights componentNames string_array {} List of component names hydrocarbonFormationVolFactorTableNames string_array {} | List of formation volume factor TableFunction names from the Functions block. diff --git a/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst b/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst index 9ad3b2e0db1..a5100a7e2f7 100644 --- a/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst +++ b/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst @@ -23,7 +23,7 @@ shotIndex integer 0 Set the current shot for tem sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. timeSourceDelay real32 -1 Source time delay (1 / f0 by default) -timeSourceFrequency real32 required Central frequency for the time source +timeSourceFrequency real32 0 Central frequency for the time source LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` ========================= ============== ========== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/ElasticFirstOrderSEM_other.rst b/src/coreComponents/schema/docs/ElasticFirstOrderSEM_other.rst index 3633753cf60..8e0d4869c8a 100644 --- a/src/coreComponents/schema/docs/ElasticFirstOrderSEM_other.rst +++ b/src/coreComponents/schema/docs/ElasticFirstOrderSEM_other.rst @@ -10,10 +10,11 @@ indexSeismoTrace integer maxStableDt real64 Value of the Maximum Stable Timestep for this solver. meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. rcvElem integer_array Element containing the receivers +receiverConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point receiverRegion integer_array Region containing the receivers -sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds +sourceConstants real64_array2d Constant part of the source for the nodes listed in m_sourceNodeIds sourceElem integer_array Element containing the sources sourceIsAccessible integer_array Flag that indicates whether the source is local to this MPI rank sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point diff --git a/src/coreComponents/schema/docs/ElasticSEM.rst b/src/coreComponents/schema/docs/ElasticSEM.rst index ea3c65a1ac4..6f0646af044 100644 --- a/src/coreComponents/schema/docs/ElasticSEM.rst +++ b/src/coreComponents/schema/docs/ElasticSEM.rst @@ -25,7 +25,7 @@ sourceForce R1Tensor {0,0,0} Force of the source: 3 re sourceMoment R2SymTensor {1,1,1,0,0,0} Moment of the source: 6 real values describing a symmetric tensor in Voigt notation.The default value is { 1, 1, 1, 0, 0, 0 } (diagonal moment, corresponding to a pure explosion). targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. timeSourceDelay real32 -1 Source time delay (1 / f0 by default) -timeSourceFrequency real32 required Central frequency for the time source +timeSourceFrequency real32 0 Central frequency for the time source LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` ========================= ============== ============= ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/ElasticSEM_other.rst b/src/coreComponents/schema/docs/ElasticSEM_other.rst index e781c2f5a58..3ee1e5dc0a7 100644 --- a/src/coreComponents/schema/docs/ElasticSEM_other.rst +++ b/src/coreComponents/schema/docs/ElasticSEM_other.rst @@ -1,25 +1,28 @@ -=========================== ============================================================================================================================================================== ======================================================================= -Name Type Description -=========================== ============================================================================================================================================================== ======================================================================= -displacementXNp1AtReceivers real32_array2d Displacement value at each receiver for each timestep (x-component) -displacementYNp1AtReceivers real32_array2d Displacement value at each receiver for each timestep (y-component) -displacementZNp1AtReceivers real32_array2d Displacement value at each receiver for each timestep (z-component) -indexSeismoTrace integer Count for output pressure at receivers -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank -receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point -sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds -sourceIsAccessible integer_array Flag that indicates whether the source is accessible to this MPI rank -sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point -sourceValue real32_array2d Source Value of the sources -useDAS integer Flag to indicate if DAS type of data will be modeled -usePML integer Flag to apply PML -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -SolverStatistics node :ref:`DATASTRUCTURE_SolverStatistics` -=========================== ============================================================================================================================================================== ======================================================================= +=========================== ============================================================================================================================================================== ================================================================================== +Name Type Description +=========================== ============================================================================================================================================================== ================================================================================== +displacementXNp1AtReceivers real32_array2d Displacement value at each receiver for each timestep (x-component) +displacementYNp1AtReceivers real32_array2d Displacement value at each receiver for each timestep (y-component) +displacementZNp1AtReceivers real32_array2d Displacement value at each receiver for each timestep (z-component) +indexSeismoTrace integer Count for output pressure at receivers +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +rcvElem integer_array Element containing the receivers +receiverConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds +receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank +receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point +receiverRegion integer_array Region containing the receivers +sourceConstants real64_array2d Constant part of the source for the nodes listed in m_sourceNodeIds in z-direction +sourceIsAccessible integer_array Flag that indicates whether the source is local to this MPI rank +sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point +sourceValue real32_array2d Source Value of the sources +useDAS integer Flag to indicate if DAS type of data will be modeled +usePML integer Flag to apply PML +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +SolverStatistics node :ref:`DATASTRUCTURE_SolverStatistics` +=========================== ============================================================================================================================================================== ================================================================================== diff --git a/src/coreComponents/schema/docs/NonlinearSolverParameters.rst b/src/coreComponents/schema/docs/NonlinearSolverParameters.rst index 47b8abb724b..4067911d031 100644 --- a/src/coreComponents/schema/docs/NonlinearSolverParameters.rst +++ b/src/coreComponents/schema/docs/NonlinearSolverParameters.rst @@ -21,6 +21,7 @@ maxAllowedResidualNorm real64 maxNumConfigurationAttempts integer 10 Max number of times that the configuration can be changed maxSubSteps integer 10 Maximum number of time sub-steps allowed for the solver maxTimeStepCuts integer 2 Max number of time step cuts +minNormalizer real64 1e-12 Value used to make sure that residual normalizers are not too small when computing residual norm. newtonMaxIter integer 5 Maximum number of iterations that are allowed in a Newton loop. newtonMinIter integer 1 Minimum number of iterations that are required before exiting the Newton loop. newtonTol real64 1e-06 The required tolerance in order to exit the Newton iteration loop. diff --git a/src/coreComponents/schema/docs/ReactiveBrine.rst b/src/coreComponents/schema/docs/ReactiveBrine.rst index d4480c90f2d..58df878c891 100644 --- a/src/coreComponents/schema/docs/ReactiveBrine.rst +++ b/src/coreComponents/schema/docs/ReactiveBrine.rst @@ -1,13 +1,14 @@ -==================== ============ ======== ============================================================================== -Name Type Default Description -==================== ============ ======== ============================================================================== -componentMolarWeight real64_array {0} Component molar weights -componentNames string_array {} List of component names -name string required A name is required for any non-unique nodes -phaseNames string_array {} List of fluid phases -phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models -==================== ============ ======== ============================================================================== +==================== ============ ======== ============================================================================================================ +Name Type Default Description +==================== ============ ======== ============================================================================================================ +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. +componentMolarWeight real64_array {0} Component molar weights +componentNames string_array {} List of component names +name string required A name is required for any non-unique nodes +phaseNames string_array {} List of fluid phases +phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models +==================== ============ ======== ============================================================================================================ diff --git a/src/coreComponents/schema/docs/ReactiveBrineThermal.rst b/src/coreComponents/schema/docs/ReactiveBrineThermal.rst index d4480c90f2d..58df878c891 100644 --- a/src/coreComponents/schema/docs/ReactiveBrineThermal.rst +++ b/src/coreComponents/schema/docs/ReactiveBrineThermal.rst @@ -1,13 +1,14 @@ -==================== ============ ======== ============================================================================== -Name Type Default Description -==================== ============ ======== ============================================================================== -componentMolarWeight real64_array {0} Component molar weights -componentNames string_array {} List of component names -name string required A name is required for any non-unique nodes -phaseNames string_array {} List of fluid phases -phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models -==================== ============ ======== ============================================================================== +==================== ============ ======== ============================================================================================================ +Name Type Default Description +==================== ============ ======== ============================================================================================================ +checkPVTTablesRanges integer 1 Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range. +componentMolarWeight real64_array {0} Component molar weights +componentNames string_array {} List of component names +name string required A name is required for any non-unique nodes +phaseNames string_array {} List of fluid phases +phasePVTParaFiles path_array required Names of the files defining the parameters of the viscosity and density models +==================== ============ ======== ============================================================================================================ diff --git a/src/coreComponents/schema/docs/SinglePhaseFVM.rst b/src/coreComponents/schema/docs/SinglePhaseFVM.rst index bd098a42ba1..c5d8ca858dd 100644 --- a/src/coreComponents/schema/docs/SinglePhaseFVM.rst +++ b/src/coreComponents/schema/docs/SinglePhaseFVM.rst @@ -3,11 +3,13 @@ ========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== Name Type Default Description ========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== +allowNegativePressure integer 1 Flag indicating if negative pressure is allowed cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] discretization string required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. isThermal integer 0 Flag indicating whether the problem is thermal or not. logLevel integer 0 Log level +maxPressureChange real64 -1 Maximum (absolute) pressure change in a Newton iteration name string required A name is required for any non-unique nodes targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. temperature real64 0 Temperature diff --git a/src/coreComponents/schema/docs/SinglePhaseHybridFVM.rst b/src/coreComponents/schema/docs/SinglePhaseHybridFVM.rst index bd098a42ba1..c5d8ca858dd 100644 --- a/src/coreComponents/schema/docs/SinglePhaseHybridFVM.rst +++ b/src/coreComponents/schema/docs/SinglePhaseHybridFVM.rst @@ -3,11 +3,13 @@ ========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== Name Type Default Description ========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== +allowNegativePressure integer 1 Flag indicating if negative pressure is allowed cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] discretization string required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. isThermal integer 0 Flag indicating whether the problem is thermal or not. logLevel integer 0 Log level +maxPressureChange real64 -1 Maximum (absolute) pressure change in a Newton iteration name string required A name is required for any non-unique nodes targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. temperature real64 0 Temperature diff --git a/src/coreComponents/schema/docs/SinglePhaseProppantFVM.rst b/src/coreComponents/schema/docs/SinglePhaseProppantFVM.rst index bd098a42ba1..c5d8ca858dd 100644 --- a/src/coreComponents/schema/docs/SinglePhaseProppantFVM.rst +++ b/src/coreComponents/schema/docs/SinglePhaseProppantFVM.rst @@ -3,11 +3,13 @@ ========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== Name Type Default Description ========================= ============ ======== ======================================================================================================================================================================================================================================================================================================================== +allowNegativePressure integer 1 Flag indicating if negative pressure is allowed cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] discretization string required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. isThermal integer 0 Flag indicating whether the problem is thermal or not. logLevel integer 0 Log level +maxPressureChange real64 -1 Maximum (absolute) pressure change in a Newton iteration name string required A name is required for any non-unique nodes targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. temperature real64 0 Temperature diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index dd35acc0008..198eb1e24cf 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1745,6 +1745,8 @@ the relative residual norm satisfies: + + @@ -2046,7 +2048,7 @@ the relative residual norm satisfies: - + @@ -2094,7 +2096,7 @@ the relative residual norm satisfies: - + @@ -2142,7 +2144,7 @@ the relative residual norm satisfies: - + @@ -2169,6 +2171,8 @@ the relative residual norm satisfies: + + @@ -2187,6 +2191,10 @@ the relative residual norm satisfies: + + + + @@ -2223,6 +2231,8 @@ the relative residual norm satisfies: + + @@ -2237,6 +2247,10 @@ the relative residual norm satisfies: + + + + @@ -2280,7 +2294,7 @@ the relative residual norm satisfies: - + @@ -2393,7 +2407,7 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> - + @@ -2445,7 +2459,7 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> - + @@ -2761,6 +2775,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -2771,6 +2787,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -2783,6 +2801,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -2793,6 +2813,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -2893,6 +2915,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -2903,6 +2927,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -3400,6 +3426,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -3476,6 +3504,8 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + @@ -3490,6 +3520,8 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + @@ -3504,6 +3536,8 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + @@ -3518,6 +3552,8 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + @@ -3566,6 +3602,8 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + @@ -3823,6 +3861,8 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + @@ -4363,6 +4403,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + @@ -4375,6 +4417,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index c7c2838827f..843c1c1b2e7 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -531,13 +531,15 @@ + + - + @@ -592,11 +594,17 @@ + + + + - + + + @@ -623,11 +631,17 @@ + + + + - + + + @@ -718,13 +732,15 @@ + + - + @@ -759,13 +775,19 @@ + + + + - + + + - + diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp index 2f93fa26a2b..44fcc1dce8a 100644 --- a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp +++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp @@ -207,13 +207,13 @@ TEST_F( AcousticWaveEquationSEMTest, SeismoTrace ) pReceivers.move( hostMemorySpace, false ); // check number of seismos and trace length - ASSERT_EQ( pReceivers.size( 1 ), 9 ); + ASSERT_EQ( pReceivers.size( 1 ), 10 ); ASSERT_EQ( pReceivers.size( 0 ), 11 ); // check seismo content. The pressure values cannot be directly checked as the problem is too small. // Since the basis is linear, check that the seismograms are nonzero (for t>0) and the seismogram at the center is equal // to the average of the others. - for( int i=0; i<11; i++ ) + for( int i = 0; i < 11; i++ ) { if( i > 0 ) { @@ -224,17 +224,17 @@ TEST_F( AcousticWaveEquationSEMTest, SeismoTrace ) { avg += pReceivers[i][r]; } - avg /=8.0; + avg /= 8.0; ASSERT_TRUE( std::abs( pReceivers[i][8] - avg ) < 0.00001 ); } // run adjoint solver - for( int i=0; i<10; i++ ) + for( int i = 0; i < 10; i++ ) { propagator->explicitStepBackward( time_n, dt, i, domain, false ); time_n += dt; } // check again the seismo content. - for( int i=0; i<11; i++ ) + for( int i = 0; i < 11; i++ ) { if( i > 0 ) { @@ -245,7 +245,7 @@ TEST_F( AcousticWaveEquationSEMTest, SeismoTrace ) { avg += pReceivers[i][r]; } - avg /=8.0; + avg /= 8.0; ASSERT_TRUE( std::abs( pReceivers[i][8] - avg ) < 0.00001 ); } } diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAcousticFirstOrder.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAcousticFirstOrder.cpp index f4167afc548..603a8bbb461 100644 --- a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAcousticFirstOrder.cpp +++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAcousticFirstOrder.cpp @@ -204,22 +204,20 @@ TEST_F( AcousticFirstOrderWaveEquationSEMTest, SeismoTrace ) uyReceivers.move( LvArray::MemorySpace::host, false ); uzReceivers.move( LvArray::MemorySpace::host, false ); - - // check number of seismos and trace length - ASSERT_EQ( pReceivers.size( 1 ), 9 ); + ASSERT_EQ( pReceivers.size( 1 ), 10 ); ASSERT_EQ( pReceivers.size( 0 ), 21 ); - ASSERT_EQ( uxReceivers.size( 1 ), 9 ); + ASSERT_EQ( uxReceivers.size( 1 ), 10 ); ASSERT_EQ( uxReceivers.size( 0 ), 21 ); - ASSERT_EQ( uyReceivers.size( 1 ), 9 ); + ASSERT_EQ( uyReceivers.size( 1 ), 10 ); ASSERT_EQ( uyReceivers.size( 0 ), 21 ); - ASSERT_EQ( uzReceivers.size( 1 ), 9 ); + ASSERT_EQ( uzReceivers.size( 1 ), 10 ); ASSERT_EQ( uzReceivers.size( 0 ), 21 ); // check seismo content. The pressure and velocity values cannot be directly checked as the problem is too small. // Since the basis is linear, check that the seismograms are nonzero (for t>0) and the seismogram at the center is equal // to the average of the others. - for( int i=0; i<21; i++ ) + for( int i = 0; i < 21; i++ ) { if( i > 0 ) { @@ -239,10 +237,10 @@ TEST_F( AcousticFirstOrderWaveEquationSEMTest, SeismoTrace ) avgUy += uyReceivers[i][r]; avgUz += uzReceivers[i][r]; } - avgP /=8.0; - avgUx /=8.0; - avgUy /=8.0; - avgUz /=8.0; + avgP /= 8.0; + avgUx /= 8.0; + avgUy /= 8.0; + avgUz /= 8.0; ASSERT_TRUE( std::abs( pReceivers[i][8] - avgP ) < 0.00001 ); ASSERT_TRUE( std::abs( uxReceivers[i][8] - avgUx ) < 0.00001 ); ASSERT_TRUE( std::abs( uyReceivers[i][8] - avgUy ) < 0.00001 );