From 639e354eeba7290ba64a79e5a0c932b884445d49 Mon Sep 17 00:00:00 2001 From: tbeltzun <129868353+tbeltzun@users.noreply.github.com> Date: Sun, 24 Sep 2023 11:07:02 +0200 Subject: [PATCH 1/4] 1st order solvers damping kernel fix --- .../AcousticFirstOrderWaveEquationSEM.cpp | 6 +- ...cousticFirstOrderWaveEquationSEMKernel.hpp | 46 +++++++------ .../ElasticFirstOrderWaveEquationSEM.cpp | 6 +- ...ElasticFirstOrderWaveEquationSEMKernel.hpp | 65 ++++++++----------- .../ElasticWaveEquationSEMKernel.hpp | 4 +- 5 files changed, 58 insertions(+), 69 deletions(-) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp index dad568c8dde..4a3c9b7aa65 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp @@ -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 >(); @@ -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, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp index 712e909fdb5..074eab3ecbc 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp @@ -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 @@ -279,42 +279,40 @@ 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 ) + localIndex const f = elemsToFaces( e, i ); + // face on the domain boundary and not on free surface + if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { - k = facesToElems( f, 1 ); - } - - real32 const alpha = 1.0 / velocity[k]; + real32 const alpha = 1.0 / velocity[e]; - constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; + constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - real64 xLocal[ numNodesPerFace ][ 3 ]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - for( localIndex i = 0; i < 3; ++i ) + 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 ); + 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 diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp index e423f910c82..4efae67e92b 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp @@ -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 >(); @@ -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, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp index c4564bbf791..5f492d5daca 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp @@ -262,8 +262,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 @@ -273,8 +273,8 @@ 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, @@ -286,46 +286,37 @@ struct DampingMatrixKernel arrayView1d< real32 > const dampingy, arrayView1d< real32 > const dampingz ) { - 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 ) + localIndex const f = elemsToFaces( e, i ); + // face on the domain boundary and not on free surface + if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { - k = facesToElems( f, 1 ); - } - - constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; + constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - real64 xLocal[ numNodesPerFace ][ 3 ]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - for( localIndex i = 0; i < 3; ++i ) + 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 localIncrementx = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][0] ) + velocityVs[k]*sqrt( faceNormal[f][1]*faceNormal[f][1] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementy = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][1] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementz = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][2] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][1]*faceNormal[f][1] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes[f][q]], localIncrementx ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes[f][q]], localIncrementy ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes[f][q]], localIncrementz ); + real32 const nx = faceNormal( f, 0 ), ny = faceNormal( f, 1 ), nz = faceNormal( f, 2 ); + for( localIndex q = 0; q < numNodesPerFace; ++q ) + { + real32 const aux = density[e] * m_finiteElement.computeDampingTerm( q, xLocal ); + real32 const localIncrementx = density[e] * (velocityVp[e] * abs( nx ) + velocityVs[e] * sqrt( pow( ny, 2 ) + pow( nz, 2 ) ) ) * aux; + real32 const localIncrementy = density[e] * (velocityVp[e] * abs( ny ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( nz, 2 ) ) ) * aux; + real32 const localIncrementz = density[e] * (velocityVp[e] * abs( nz ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( ny, 2 ) ) ) * aux; + + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes( f, q )], localIncrementx ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes( f, q )], localIncrementy ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes( f, q )], localIncrementz ); + } } } } ); // end loop over element diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp index a776dec8a6c..bfe020d7005 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp @@ -286,8 +286,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 From 1b6444dd9e96e300a0f9ebec8772824c0d8b9a58 Mon Sep 17 00:00:00 2001 From: tbeltzun <129868353+tbeltzun@users.noreply.github.com> Date: Sun, 10 Sep 2023 21:33:34 +0200 Subject: [PATCH 2/4] `DampingMatrixKernel` faces fix for 2nd order equations --- .../AcousticWaveEquationSEM.cpp | 31 ++++---- .../AcousticWaveEquationSEMKernel.hpp | 56 +++++++-------- .../ElasticWaveEquationSEM.cpp | 50 ++++++------- .../ElasticWaveEquationSEMKernel.hpp | 72 ++++++++----------- 4 files changed, 95 insertions(+), 114 deletions(-) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp index f8fef564ce3..97bab9b5eec 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp @@ -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(); @@ -283,28 +284,29 @@ 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 >(); + arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensityA >(); /// 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, @@ -312,10 +314,9 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() { 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, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp index c0a9044f7b6..006024274c0 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp @@ -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 } @@ -252,7 +252,7 @@ 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] nodeCoords coordinates of the nodes * @param[in] facesToElems map from faces to elements * @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 @@ -264,8 +264,8 @@ 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, @@ -273,37 +273,33 @@ struct DampingMatrixKernel 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 diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp index ca71b4c7989..bce5887a4e1 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp @@ -525,14 +525,9 @@ void ElasticWaveEquationSEM::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 X = nodeManager.getField< fields::referencePosition32 >().toViewConst(); - arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - - /// get face to node map - ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst(); // mass matrix to be computed in this function arrayView1d< real32 > const mass = nodeManager.getField< fields::MassVector >(); @@ -546,40 +541,39 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() dampingz.zero(); /// get array of indicators: 1 if face is on the free surface; 0 otherwise - arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >(); + arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >(); + arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator(); + ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); + arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - 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(); - arrayView1d< real32 > const density = elementSubRegion.getField< fields::MediumDensity >(); - arrayView1d< real32 > const velocityVp = elementSubRegion.getField< fields::MediumVelocityVp >(); - arrayView1d< real32 > const velocityVs = elementSubRegion.getField< fields::MediumVelocityVs >(); + arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); - finiteElement::FiniteElementBase const & - fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensityE >(); + arrayView1d< real32 const > const velocityVp = elementSubRegion.getField< fields::MediumVelocityVp >(); + arrayView1d< real32 const > const velocityVs = elementSubRegion.getField< fields::MediumVelocityVs >(); - finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, - [&] - ( auto const finiteElement ) + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { using FE_TYPE = TYPEOFREF( finiteElement ); elasticWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement ); - kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), - X, + nodeCoords, elemsToNodes, density, mass ); elasticWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement ); - - kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(), - X, - facesToElements, + kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), + nodeCoords, + elemsToFaces, facesToNodes, facesDomainBoundaryIndicator, freeSurfaceFaceIndicator, @@ -608,9 +602,9 @@ void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartit arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >(); arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >(); arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >(); - arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >(); - arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >(); - arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >(); + arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >(); + arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >(); + arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >(); arrayView1d< real32 > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >(); arrayView1d< real32 > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >(); arrayView1d< real32 > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >(); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp index bfe020d7005..5f9ad660db3 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp @@ -245,7 +245,7 @@ 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; @@ -256,14 +256,14 @@ struct MassMatrixKernel { 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 = density[k] * m_finiteElement.computeMassTerm( q, xLocal ); - RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement ); + real32 const localIncrement = density[e] * m_finiteElement.computeMassTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( e, q )], localIncrement ); } } ); // end loop over element } @@ -298,8 +298,8 @@ struct DampingMatrixKernel //std::enable_if_t< geos::is_sem_formulation< std::remove_cv_t< FE_TYPE_ > >::value, void > 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, @@ -311,49 +311,39 @@ struct DampingMatrixKernel arrayView1d< real32 > const dampingy, arrayView1d< real32 > const dampingz ) { - 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 ); - } - - 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 localIncrementx = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][0] ) + velocityVs[k]*sqrt( faceNormal[f][1]*faceNormal[f][1] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementy = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][1] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementz = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][2] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][1]*faceNormal[f][1] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes[f][q]], localIncrementx ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes[f][q]], localIncrementy ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes[f][q]], localIncrementz ); + real32 const nx = faceNormal( f, 0 ), ny = faceNormal( f, 1 ), nz = faceNormal( f, 2 ); + for( localIndex q = 0; q < numNodesPerFace; ++q ) + { + real32 const aux = density[e] * m_finiteElement.computeDampingTerm( q, xLocal ); + real32 const localIncrementx = aux * ( velocityVp[e] * abs( nx ) + velocityVs[e] * sqrt( pow( ny, 2 ) + pow( nz, 2 ) ) ); + real32 const localIncrementy = aux * ( velocityVp[e] * abs( ny ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( nz, 2 ) ) ); + real32 const localIncrementz = aux * ( velocityVp[e] * abs( nz ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( ny, 2 ) ) ); + + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes( f, q )], localIncrementx ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes( f, q )], localIncrementy ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes( f, q )], localIncrementz ); + } } } - } ); // end loop over element + } ); } /// The finite element space/discretization object for the element type in the subRegion From 42853a65ea1160a6a3d25c218263f707dce64fb3 Mon Sep 17 00:00:00 2001 From: tbeltzun <129868353+tbeltzun@users.noreply.github.com> Date: Fri, 29 Sep 2023 16:24:09 +0200 Subject: [PATCH 3/4] fix docstrings - format --- .../AcousticFirstOrderWaveEquationSEMKernel.hpp | 6 ++---- .../wavePropagation/AcousticWaveEquationSEM.cpp | 2 +- .../wavePropagation/AcousticWaveEquationSEMKernel.hpp | 2 +- .../ElasticFirstOrderWaveEquationSEMKernel.hpp | 3 +-- 4 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp index 074eab3ecbc..cc01a3ba96a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp @@ -295,10 +295,7 @@ struct DampingMatrixKernel // face on the domain boundary and not on free surface if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { - real32 const alpha = 1.0 / velocity[e]; - constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - real64 xLocal[ numNodesPerFace ][ 3 ]; for( localIndex a = 0; a < numNodesPerFace; ++a ) { @@ -308,6 +305,7 @@ struct DampingMatrixKernel } } + real32 const alpha = 1.0 / velocity[e]; for( localIndex q = 0; q < numNodesPerFace; ++q ) { real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal ); @@ -315,7 +313,7 @@ struct DampingMatrixKernel } } } - } ); // end loop over element + } ); } /// The finite element space/discretization object for the element type in the subRegion diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp index 97bab9b5eec..14d1efb9750 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp @@ -299,7 +299,7 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() /// Partial gradient if gradient as to be computed arrayView1d< real32 > grad = elementSubRegion.getField< fields::PartialGradient >(); grad.zero(); - + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { using FE_TYPE = TYPEOFREF( finiteElement ); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp index 006024274c0..b9815ce1b94 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp @@ -253,7 +253,7 @@ struct DampingMatrixKernel * @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] facesToElems map from faces to elements + * @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 diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp index 5f492d5daca..973ca5c8403 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp @@ -295,7 +295,6 @@ struct DampingMatrixKernel if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - real64 xLocal[ numNodesPerFace ][ 3 ]; for( localIndex a = 0; a < numNodesPerFace; ++a ) { @@ -319,7 +318,7 @@ struct DampingMatrixKernel } } } - } ); // end loop over element + } ); } /// The finite element space/discretization object for the element type in the subRegion From aa473161236036340583d3af1663b573e59400d8 Mon Sep 17 00:00:00 2001 From: tbeltzun <129868353+tbeltzun@users.noreply.github.com> Date: Fri, 29 Sep 2023 16:31:14 +0200 Subject: [PATCH 4/4] fix unrelated change --- .../physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp | 2 +- .../physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp index 14d1efb9750..6edd7133e4d 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp @@ -294,7 +294,7 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView1d< real32 const > const velocity = elementSubRegion.getField< fields::MediumVelocity >(); - arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensityA >(); + arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensity >(); /// Partial gradient if gradient as to be computed arrayView1d< real32 > grad = elementSubRegion.getField< fields::PartialGradient >(); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp index bce5887a4e1..53ff49f3d42 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp @@ -555,7 +555,7 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); - arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensityE >(); + arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensity >(); arrayView1d< real32 const > const velocityVp = elementSubRegion.getField< fields::MediumVelocityVp >(); arrayView1d< real32 const > const velocityVs = elementSubRegion.getField< fields::MediumVelocityVs >();