Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rewrite wave propagation damping kernels #2730

Merged
merged 8 commits into from
Oct 11, 2023
Merged

rewrite wave propagation damping kernels #2730

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
639e354
1st order solvers damping kernel fix
tbeltzun Sep 24, 2023
1b6444d
`DampingMatrixKernel` faces fix for 2nd order equations
tbeltzun Sep 10, 2023
42853a6
fix docstrings - format
tbeltzun Sep 29, 2023
aa47316
fix unrelated change
tbeltzun Sep 29, 2023
54d1c63
Merge branch 'develop' into bugfix/tbeltzun/damping-faces
tbeltzun Oct 9, 2023
0f8ccd9
Merge branch 'develop' into bugfix/tbeltzun/damping-faces
sframba Oct 9, 2023
0b99acd
Merge branch 'develop' into bugfix/tbeltzun/damping-faces
sframba Oct 9, 2023
675483a
Merge branch 'develop' into bugfix/tbeltzun/damping-faces
tbeltzun Oct 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -323,7 +323,7 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
CellElementSubRegion & elementSubRegion )
{
arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
arrayView2d< localIndex const > const facesToElements = faceManager.elementList();
arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList();
arrayView1d< real32 const > const velocity = elementSubRegion.getField< wavesolverfields::MediumVelocity >();
arrayView1d< real32 const > const density = elementSubRegion.getField< wavesolverfields::MediumDensity >();

Expand All @@ -344,9 +344,9 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro

acousticFirstOrderWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement );

kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
X,
facesToElements,
elemsToFaces,
facesToNodes,
facesDomainBoundaryIndicator,
freeSurfaceFaceIndicator,
Expand Down
50 changes: 23 additions & 27 deletions ...coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,8 +268,8 @@ struct DampingMatrixKernel
* @tparam EXEC_POLICY the execution policy
* @tparam ATOMIC_POLICY the atomic policy
* @param[in] size the number of cells in the subRegion
* @param[in] X coordinates of the nodes
* @param[in] facesToElems map from faces to elements
* @param[in] nodeCoords coordinates of the nodes
* @param[in] elemsToFaces map from elements to faces
* @param[in] facesToNodes map from face to nodes
* @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise
* @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise
Expand All @@ -279,45 +279,41 @@ struct DampingMatrixKernel
template< typename EXEC_POLICY, typename ATOMIC_POLICY >
void
launch( localIndex const size,
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
arrayView2d< localIndex const > const facesToElems,
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
arrayView2d< localIndex const > const elemsToFaces,
ArrayOfArraysView< localIndex const > const facesToNodes,
arrayView1d< integer const > const facesDomainBoundaryIndicator,
arrayView1d< localIndex const > const freeSurfaceFaceIndicator,
arrayView1d< real32 const > const velocity,
arrayView1d< real32 > const damping )
{
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f )
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
// face on the domain boundary and not on free surface
if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i )
{
localIndex k = facesToElems( f, 0 );
if( k < 0 )
{
k = facesToElems( f, 1 );
}

real32 const alpha = 1.0 / velocity[k];

constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;

real64 xLocal[ numNodesPerFace ][ 3 ];
for( localIndex a = 0; a < numNodesPerFace; ++a )
localIndex const f = elemsToFaces( e, i );
// face on the domain boundary and not on free surface
if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
{
for( localIndex i = 0; i < 3; ++i )
constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
real64 xLocal[ numNodesPerFace ][ 3 ];
for( localIndex a = 0; a < numNodesPerFace; ++a )
{
xLocal[a][i] = X( facesToNodes( f, a ), i );
for( localIndex d = 0; d < 3; ++d )
{
xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d );
}
}
}

for( localIndex q = 0; q < numNodesPerFace; ++q )
{
real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes[f][q]], localIncrement );
real32 const alpha = 1.0 / velocity[e];
for( localIndex q = 0; q < numNodesPerFace; ++q )
{
real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes( f, q )], localIncrement );
}
}
}
} ); // end loop over element
} );
}

/// The finite element space/discretization object for the element type in the subRegion
Expand Down
29 changes: 15 additions & 14 deletions src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -260,10 +260,11 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()

NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
ElementRegionManager & elemManager = mesh.getElemManager();

/// get the array of indicators: 1 if the face is on the boundary; 0 otherwise
arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator();
arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X32 = nodeManager.getField< fields::referencePosition32 >().toViewConst();
arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst();

/// get face to node map
ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();
Expand All @@ -283,39 +284,39 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
/// get array of indicators: 1 if face is on the free surface; 0 otherwise
arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >();

mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
CellElementSubRegion & elementSubRegion )
elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
CellElementSubRegion & elementSubRegion )
{
finiteElement::FiniteElementBase const &
fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );

arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
arrayView2d< localIndex const > const facesToElements = faceManager.elementList();
arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList();

arrayView1d< real32 const > const velocity = elementSubRegion.getField< fields::MediumVelocity >();
arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensity >();

/// Partial gradient if gradient as to be computed
arrayView1d< real32 > grad = elementSubRegion.getField< fields::PartialGradient >();
{
grad.zero();
}
finiteElement::FiniteElementBase const &
fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
grad.zero();

finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );

acousticWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement );
kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
X32,
nodeCoords,
elemsToNodes,
velocity,
density,
mass );
{
GEOS_MARK_SCOPE( DampingMatrixKernel );
acousticWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement );

kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
X32,
facesToElements,
kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
nodeCoords,
elemsToFaces,
facesToNodes,
facesDomainBoundaryIndicator,
freeSurfaceFaceIndicator,
Expand Down
58 changes: 27 additions & 31 deletions src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,25 +211,25 @@ struct MassMatrixKernel
arrayView1d< real32 > const mass )

{
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints;

real32 const invC2 = 1.0 / ( density[k] * velocity[k] * velocity[k] );
real32 const invC2 = 1.0 / ( density[e] * velocity[e] * velocity[e] );
real64 xLocal[ numNodesPerElem ][ 3 ];
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
for( localIndex i = 0; i < 3; ++i )
{
xLocal[a][i] = X( elemsToNodes( k, a ), i );
xLocal[a][i] = X( elemsToNodes( e, a ), i );
}
}

for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q )
{
real32 const localIncrement = invC2 * m_finiteElement.computeMassTerm( q, xLocal );
RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement );
RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( e, q )], localIncrement );
}
} ); // end loop over element
}
Expand All @@ -252,8 +252,8 @@ struct DampingMatrixKernel
* @tparam EXEC_POLICY the execution policy
* @tparam ATOMIC_POLICY the atomic policy
* @param[in] size the number of cells in the subRegion
* @param[in] X coordinates of the nodes
* @param[in] facesToElems map from faces to elements
* @param[in] nodeCoords coordinates of the nodes
* @param[in] elemsToFaces map from elements to faces
* @param[in] facesToNodes map from face to nodes
* @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise
* @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise
Expand All @@ -264,46 +264,42 @@ struct DampingMatrixKernel
template< typename EXEC_POLICY, typename ATOMIC_POLICY >
void
launch( localIndex const size,
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
arrayView2d< localIndex const > const facesToElems,
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
arrayView2d< localIndex const > const elemsToFaces,
ArrayOfArraysView< localIndex const > const facesToNodes,
arrayView1d< integer const > const facesDomainBoundaryIndicator,
arrayView1d< localIndex const > const freeSurfaceFaceIndicator,
arrayView1d< real32 const > const velocity,
arrayView1d< real32 const > const density,
arrayView1d< real32 > const damping )
{
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f )
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
// face on the domain boundary and not on free surface
if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i )
{
localIndex k = facesToElems( f, 0 );
if( k < 0 )
{
k = facesToElems( f, 1 );
}

real32 const alpha = 1.0 / (density[k] * velocity[k]);

constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;

real64 xLocal[ numNodesPerFace ][ 3 ];
for( localIndex a = 0; a < numNodesPerFace; ++a )
localIndex const f = elemsToFaces( e, i );
// face on the domain boundary and not on free surface
if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
{
for( localIndex i = 0; i < 3; ++i )
constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
real64 xLocal[ numNodesPerFace ][ 3 ];
for( localIndex a = 0; a < numNodesPerFace; ++a )
{
xLocal[a][i] = X( facesToNodes( f, a ), i );
for( localIndex d = 0; d < 3; ++d )
{
xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d );
}
}
}
real32 const alpha = 1.0 / (density[e] * velocity[e]);

for( localIndex q = 0; q < numNodesPerFace; ++q )
{
real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes[f][q]], localIncrement );
for( localIndex q = 0; q < numNodesPerFace; ++q )
{
real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes( f, q )], localIncrement );
}
}
}
} ); // end loop over element
} );
}

/// The finite element space/discretization object for the element type in the subRegion
Expand Down
6 changes: 3 additions & 3 deletions src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -383,7 +383,7 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou
{

arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
arrayView2d< localIndex const > const facesToElements = faceManager.elementList();
arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList();
arrayView1d< real32 > const density = elementSubRegion.getField< wavesolverfields::MediumDensity >();
arrayView1d< real32 > const velocityVp = elementSubRegion.getField< wavesolverfields::MediumVelocityVp >();
arrayView1d< real32 > const velocityVs = elementSubRegion.getField< wavesolverfields::MediumVelocityVs >();
Expand All @@ -407,9 +407,9 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou

elasticFirstOrderWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement );

kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
X,
facesToElements,
elemsToFaces,
facesToNodes,
facesDomainBoundaryIndicator,
freeSurfaceFaceIndicator,
Expand Down
Loading