Skip to content

Commit

Permalink
rewrite seismo traces - rework WaveSolverBase structure (#2729)
Browse files Browse the repository at this point in the history 
* rewrite seismo traces - rework `WaveSolverBase` structure
* format
* fix tests
* update submodules
* update PR for VTI
* consistency
* update submodules
* updated integratedTests submodule
* Updating schema
* Updating the integratedTests hash

---------

Co-authored-by: Julien Rene Thomas Besset <[email protected]>
Co-authored-by: acitrain <[email protected]>
Co-authored-by: Pavel Tomin <[email protected]>
Co-authored-by: Nicola Castelletto <[email protected]>
  • Loading branch information
@tbeltzun @acitrain
@paveltomin @castelletto1
5 people authored Nov 9, 2023
1 parent 23abbb8 commit 625714b
Show file tree
Hide file tree
Showing 45 changed files with 596 additions and 728 deletions.
2 changes: 1 addition & 1 deletion integratedTests
Submodule integratedTests updated 265 files
104 changes: 27 additions & 77 deletions src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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 );
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -354,6 +342,8 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
} );
} );
} );

WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
}


Expand Down Expand Up @@ -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 );
Expand All @@ -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;
}

Expand All @@ -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" );
Expand Down
42 changes: 1 addition & 41 deletions src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -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
Expand All @@ -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;

Expand Down Expand Up @@ -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;

};


Expand Down
52 changes: 12 additions & 40 deletions src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 &,
Expand Down Expand Up @@ -315,6 +314,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
} );
} );

WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
}

void AcousticVTIWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain )
Expand Down Expand Up @@ -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 )
{
Expand Down Expand Up @@ -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 */
15 changes: 0 additions & 15 deletions src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*/
Expand Down
Loading

0 comments on commit 625714b

Please sign in to comment.