Skip to content

Commit

Permalink
feat: WaveEquation : new acoustic gradient formulation (#3361)
Browse files Browse the repository at this point in the history
  • Loading branch information
@rmadec-cs
rmadec-cs authored Dec 4, 2024
1 parent dc6a6be commit b56762c
Show file tree
Hide file tree
Showing 8 changed files with 569 additions and 50 deletions.
2 changes: 1 addition & 1 deletion .integrated_tests.yaml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
baselines:
bucket: geosx
baseline: integratedTests/baseline_integratedTests-pr3393-9089-c187f5d
baseline: integratedTests/baseline_integratedTests-pr3361-9139-2fc4131
allow_fail:
all: ''
streak: ''
6 changes: 5 additions & 1 deletion BASELINE_NOTES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,11 @@ This file is designed to track changes to the integrated test baselines.
Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).

PR #3393 (2024-12-2)
PR #3361 (2024-12-03)
=====================
Baseline diffs after reimplementation of wave equation acoustic gradient for velocity and density parameters: new field "partialGradient2" and "pressureForward" field replacing "pressureDoubleDerivative".

PR #3393 (2024-12-02)
=====================
Fix netToGross bug.

Expand Down
92 changes: 58 additions & 34 deletions ...Solvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
nodeManager.registerField< acousticfields::Pressure_nm1,
acousticfields::Pressure_n,
acousticfields::Pressure_np1,
acousticfields::PressureDoubleDerivative,
acousticfields::PressureForward,
acousticfields::ForcingRHS,
acousticfields::AcousticMassVector,
acousticfields::DampingVector,
Expand Down Expand Up @@ -106,6 +106,7 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
subRegion.registerField< acousticfields::AcousticVelocity >( getName() );
subRegion.registerField< acousticfields::AcousticDensity >( getName() );
subRegion.registerField< acousticfields::PartialGradient >( getName() );
subRegion.registerField< acousticfields::PartialGradient2 >( getName() );
} );

} );
Expand Down Expand Up @@ -298,6 +299,8 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
arrayView1d< real32 > grad = elementSubRegion.getField< acousticfields::PartialGradient >();

grad.zero();
arrayView1d< real32 > grad2 = elementSubRegion.getField< acousticfields::PartialGradient2 >();
grad2.zero();

finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
Expand Down Expand Up @@ -881,43 +884,35 @@ real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n,
{
NodeManager & nodeManager = mesh.getNodeManager();

arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticfields::Pressure_nm1 >();
arrayView1d< real32 > const p_n = nodeManager.getField< acousticfields::Pressure_n >();
arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticfields::Pressure_np1 >();

if( computeGradient && cycleNumber >= 0 )
{

arrayView1d< real32 > const p_dt2 = nodeManager.getField< acousticfields::PressureDoubleDerivative >();

if( m_enableLifo )
{
if( !m_lifo )
{
int const rank = MpiWrapper::commRank( MPI_COMM_GEOS );
std::string lifoPrefix = GEOS_FMT( "lifo/rank_{:05}/pdt2_shot{:06}", rank, m_shotIndex );
m_lifo = std::make_unique< LifoStorage< real32, localIndex > >( lifoPrefix, p_dt2, m_lifoOnDevice, m_lifoOnHost, m_lifoSize );
std::string lifoPrefix = GEOS_FMT( "lifo/rank_{:05}/p_forward_shot{:06}", rank, m_shotIndex );
m_lifo = std::make_unique< LifoStorage< real32, localIndex > >( lifoPrefix, p_n, m_lifoOnDevice, m_lifoOnHost, m_lifoSize );
}

m_lifo->pushWait();
}
forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const nodeIdx )
{
p_dt2[nodeIdx] = (p_np1[nodeIdx] - 2*p_n[nodeIdx] + p_nm1[nodeIdx]) / pow( dt, 2 );
} );

if( m_enableLifo )
{
// Need to tell LvArray data is on GPU to avoir HtoD copy
p_dt2.move( LvArray::MemorySpace::cuda, false );
m_lifo->pushAsync( p_dt2 );
p_n.move( LvArray::MemorySpace::cuda, false );
m_lifo->pushAsync( p_n );
}
else
{
GEOS_MARK_SCOPE ( DirectWrite );
p_dt2.move( LvArray::MemorySpace::host, false );
p_n.move( LvArray::MemorySpace::host, false );
int const rank = MpiWrapper::commRank( MPI_COMM_GEOS );
std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressuredt2_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber );
std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressure_forward_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber );
int lastDirSeparator = fileName.find_last_of( "/\\" );
std::string dirName = fileName.substr( 0, lastDirSeparator );
if( string::npos != (size_t)lastDirSeparator && !directoryExists( dirName ))
Expand All @@ -929,7 +924,7 @@ real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n,
GEOS_THROW_IF( !wf,
getDataContext() << ": Could not open file "<< fileName << " for writing",
InputError );
wf.write( (char *)&p_dt2[0], p_dt2.size()*sizeof( real32 ) );
wf.write( (char *)&p_n[0], p_n.size()*sizeof( real32 ) );
wf.close( );
GEOS_THROW_IF( !wf.good(),
getDataContext() << ": An error occured while writing "<< fileName,
Expand Down Expand Up @@ -959,12 +954,18 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
{
NodeManager & nodeManager = mesh.getNodeManager();

arrayView1d< real32 const > const mass = nodeManager.getField< acousticfields::AcousticMassVector >();

arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticfields::Pressure_nm1 >();
arrayView1d< real32 > const p_n = nodeManager.getField< acousticfields::Pressure_n >();
arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticfields::Pressure_np1 >();

//// Compute q_dt2 and store it in p_nm1
SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst();
forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
{
localIndex const a = solverTargetNodesSet[n];
p_nm1[a] = (p_np1[a] - 2*p_n[a] + p_nm1[a]) / pow( dt, 2 );
} );

EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() );
int const maxCycle = int(round( maxTime / dt ));
Expand All @@ -973,11 +974,11 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
{
ElementRegionManager & elemManager = mesh.getElemManager();

arrayView1d< real32 > const p_dt2 = nodeManager.getField< acousticfields::PressureDoubleDerivative >();
arrayView1d< real32 > const p_forward = nodeManager.getField< acousticfields::PressureForward >();

if( m_enableLifo )
{
m_lifo->pop( p_dt2 );
m_lifo->pop( p_forward );
if( m_lifo->empty() )
delete m_lifo.release();
}
Expand All @@ -986,37 +987,60 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
GEOS_MARK_SCOPE ( DirectRead );

int const rank = MpiWrapper::commRank( MPI_COMM_GEOS );
std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressuredt2_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber );
std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressure_forward_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber );
std::ifstream wf( fileName, std::ios::in | std::ios::binary );
GEOS_THROW_IF( !wf,
getDataContext() << ": Could not open file "<< fileName << " for reading",
InputError );

p_dt2.move( LvArray::MemorySpace::host, true );
wf.read( (char *)&p_dt2[0], p_dt2.size()*sizeof( real32 ) );
p_forward.move( LvArray::MemorySpace::host, true );
wf.read( (char *)&p_forward[0], p_forward.size()*sizeof( real32 ) );
wf.close( );
remove( fileName.c_str() );
}
elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
CellElementSubRegion & elementSubRegion )
{
arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfields::AcousticVelocity >();
arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst();
arrayView1d< real32 > grad = elementSubRegion.getField< acousticfields::PartialGradient >();
arrayView1d< real32 > grad2 = elementSubRegion.getField< acousticfields::PartialGradient2 >();
arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList();
constexpr localIndex numNodesPerElem = 8;
arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank();
GEOS_MARK_SCOPE ( updatePartialGradient );
forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const eltIdx )

//COMPUTE GRADIENTS with respect to K=1/rho*c2 (grad) and b=1/rho (grad2)
finiteElement::FiniteElementBase const &
fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );

finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
if( elemGhostRank[eltIdx]<0 )
{
for( localIndex i = 0; i < numNodesPerElem; ++i )
{
localIndex nodeIdx = elemsToNodes[eltIdx][i];
grad[eltIdx] += (-2/velocity[eltIdx]) * mass[nodeIdx]/8.0 * (p_dt2[nodeIdx] * p_n[nodeIdx]);
}
}
using FE_TYPE = TYPEOFREF( finiteElement );

AcousticMatricesSEM::GradientKappaBuoyancy< FE_TYPE > kernelG( finiteElement );
kernelG.template computeGradient< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
nodeCoords,
elemsToNodes,
elemGhostRank,
p_nm1,
p_n,
p_forward,
grad,
grad2 );


} );

// // Change of variables to return grad with respect to c and rho
// arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfields::AcousticVelocity >();
// arrayView1d< real32 const > const density = elementSubRegion.getField< acousticfields::AcousticDensity >();
// forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const eltIdx )
// {
// if( elemGhostRank[eltIdx]<0 )
// {
// grad2[eltIdx] = -1/(pow(density[eltIdx]*velocity[eltIdx],2)) * grad[eltIdx] - 1/pow(density[eltIdx],2) * grad2[eltIdx];
// grad[eltIdx]= -2/(density[eltIdx]*pow(velocity[eltIdx],3)) * grad[eltIdx];
// }
// } );
} );
}

Expand Down
14 changes: 11 additions & 3 deletions src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,13 +59,13 @@ DECLARE_FIELD( Pressure_np1,
WRITE_AND_READ,
"Scalar pressure at time n+1." );

DECLARE_FIELD( PressureDoubleDerivative,
"pressureDoubleDerivative",
DECLARE_FIELD( PressureForward,
"pressureForward",
array1d< real32 >,
0,
NOPLOT,
WRITE_AND_READ,
"Double derivative of the pressure for each node to compute the gradient" );
"Pressure field from forward pass on each node to compute the gradient" );

DECLARE_FIELD( Velocity_x,
"velocity_x",
Expand Down Expand Up @@ -99,6 +99,14 @@ DECLARE_FIELD( PartialGradient,
WRITE_AND_READ,
"Partiel gradient computed during backward propagation" );

DECLARE_FIELD( PartialGradient2,
"partialGradient2",
array1d< real32 >,
0,
NOPLOT,
WRITE_AND_READ,
"Partial gradient for density/velocity computed during backward propagation" );

DECLARE_FIELD( ForcingRHS,
"rhs",
array1d< real32 >,
Expand Down
63 changes: 63 additions & 0 deletions ...mponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,70 @@ struct AcousticMatricesSEM

};

template< typename FE_TYPE >
struct GradientKappaBuoyancy
{

GradientKappaBuoyancy( FE_TYPE const & finiteElement )
: m_finiteElement( finiteElement )
{}

/**
* @brief Launch the computation of the 2 gradients relative to the coeff of the wave equation K=1/rho*c2 and b=1/rho
* @tparam EXEC_POLICY the execution policy
* @tparam ATOMIC_POLICY the atomic policy
* @param[in] size the number of cells in the subRegion
* @param[in] nodeCoords coordinates of the nodes
* @param[in] elemsToNodes map from element to nodes
* @param[in] q_dt2 second order derivative in time of backward
* @param[in] q_n current time step of backward
* @param[in] p_n current time step of forward
* @param[out] grad first part of gradient vector with respect to K=1/rho*c2
* @param[out] grad2 second part of gradient vector with respact to b=1/rho
*/
template< typename EXEC_POLICY, typename ATOMIC_POLICY >
void
computeGradient( localIndex const size,
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes,
arrayView1d< integer const > const elemGhostRank,
arrayView1d< real32 const > const q_dt2,
arrayView1d< real32 const > const q_n,
arrayView1d< real32 const > const p_n,
arrayView1d< real32 > const grad,
arrayView1d< real32 > const grad2 )

{
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
if( elemGhostRank[e]<0 )
{
// only the eight corners of the mesh cell are needed to compute the Jacobian
real64 xLocal[ 8 ][ 3 ];
for( localIndex a = 0; a < 8; ++a )
{
localIndex const nodeIndex = elemsToNodes( e, FE_TYPE::meshIndexToLinearIndex3D( a ) );
for( localIndex i = 0; i < 3; ++i )
{
xLocal[a][i] = nodeCoords( nodeIndex, i );
}
}
constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints;
for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q )
{
localIndex nodeIdx = elemsToNodes( e, q );
grad[e] += q_dt2[nodeIdx] * p_n[nodeIdx] * m_finiteElement.computeMassTerm( q, xLocal );
m_finiteElement.template computeStiffnessTerm( q, xLocal, [&] ( const int i, const int j, const real64 val )
{
grad2[e] += val* q_n[elemsToNodes( e, j )] * p_n[elemsToNodes( e, i )];
} );
}
}
} ); // end loop over element
}
/// The finite element space/discretization object for the element type in the subRegion
FE_TYPE const & m_finiteElement;
};

};

Expand Down
34 changes: 24 additions & 10 deletions src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -484,9 +484,19 @@ void WaveSolverBase::computeTargetNodeSet( arrayView2d< localIndex const, cells:

void WaveSolverBase::incrementIndexSeismoTrace( real64 const time_n )
{
while( (m_dtSeismoTrace * m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
if( m_forward )
{
while( (m_dtSeismoTrace * m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
{
m_indexSeismoTrace++;
}
}
else
{
m_indexSeismoTrace++;
while( (m_dtSeismoTrace * m_indexSeismoTrace) >= (time_n - epsilonLoc) && m_indexSeismoTrace > 0 )
{
m_indexSeismoTrace--;
}
}
}

Expand Down Expand Up @@ -515,12 +525,14 @@ void WaveSolverBase::computeAllSeismoTraces( real64 const time_n,
if( m_nsamplesSeismoTrace == 0 )
return;
integer const dir = m_forward ? +1 : -1;
for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ )
integer const beginIndex = m_forward ? m_indexSeismoTrace : m_nsamplesSeismoTrace-m_indexSeismoTrace;
for( localIndex iSeismo = beginIndex; iSeismo < m_nsamplesSeismoTrace; iSeismo++ )
{
real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo);
if( dir * timeSeismo > dir * (time_n + epsilonLoc) )
localIndex seismoIndex = m_forward ? iSeismo : m_nsamplesSeismoTrace-iSeismo;
real64 const timeSeismo = m_dtSeismoTrace * seismoIndex;
if( dir * timeSeismo > dir * time_n + epsilonLoc )
break;
WaveSolverUtils::computeSeismoTrace( time_n, dir * dt, timeSeismo, iSeismo, m_receiverNodeIds,
WaveSolverUtils::computeSeismoTrace( time_n, dir * dt, timeSeismo, seismoIndex, m_receiverNodeIds,
m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers, coeffs, add );
}
}
Expand All @@ -535,12 +547,14 @@ void WaveSolverBase::compute2dVariableAllSeismoTraces( localIndex const regionIn
if( m_nsamplesSeismoTrace == 0 )
return;
integer const dir = m_forward ? +1 : -1;
for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ )
integer const beginIndex = m_forward ? m_indexSeismoTrace : m_nsamplesSeismoTrace-m_indexSeismoTrace;
for( localIndex iSeismo = beginIndex; iSeismo < m_nsamplesSeismoTrace; iSeismo++ )
{
real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo);
if( dir * timeSeismo > dir * (time_n + epsilonLoc))
localIndex seismoIndex = m_forward ? iSeismo : m_nsamplesSeismoTrace-iSeismo;
real64 const timeSeismo = m_dtSeismoTrace * seismoIndex;
if( dir * timeSeismo > dir * time_n + epsilonLoc )
break;
WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dir * dt, regionIndex, m_receiverRegion, timeSeismo, iSeismo, m_receiverElem,
WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dir * dt, regionIndex, m_receiverRegion, timeSeismo, seismoIndex, m_receiverElem,
m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers );
}
}
Expand Down
4 changes: 3 additions & 1 deletion src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,9 @@ set( gtest_geosx_tests
testWavePropagationDAS.cpp
testWavePropagationElasticVTI.cpp
testWavePropagationAttenuation.cpp
testWavePropagationAcousticFirstOrder.cpp )
testWavePropagationAcousticFirstOrder.cpp
testWavePropagationAdjoint1.cpp
)

set( tplDependencyList ${parallelDeps} gtest )

Expand Down
Loading

0 comments on commit b56762c

Please sign in to comment.