From e43678a1b24f6a7f262252ca2e3e769a35afd3ba Mon Sep 17 00:00:00 2001 From: j0405284 Date: Wed, 21 Aug 2024 17:20:05 +0200 Subject: [PATCH 1/6] fixed a bug due to PR 3202 concerning the elastic solver at higher orders --- .../Qk_Hexahedron_Lagrange_GaussLobatto.hpp | 96 +++++++++++++++++++ .../PrecomputeSourcesAndReceiversKernel.hpp | 20 ++-- 2 files changed, 106 insertions(+), 10 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp index f7f384241be..8058ec36baa 100644 --- a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp @@ -328,6 +328,51 @@ class Qk_Hexahedron_Lagrange_GaussLobatto final : public FiniteElementBase real64 const (&X)[numNodes][3], StackVariables const & stack, real64 ( &gradN )[numNodes][3] ); + /** + * @brief Calculate the shape functions derivatives wrt the physical + * coordinates. + * @param q Index of the quadrature point. + * @param X Array containing the coordinates of the mesh corners. + * @param gradN Array to contain the shape function derivatives for all + * support points at the coordinates of the quadrature point @p q. + * @return The determinant of the parent/physical transformation matrix. + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + static real64 calcGradNWithCorners( localIndex const q, + real64 const (&X)[8][3], + real64 ( &gradN )[numNodes][3] ); + /** + * @brief Calculate the shape functions derivatives wrt the physical + * coordinates at a single point. + * @param[in] coords The parent coordinates at which to evaluate the shape function value + * @param[in] X Array containing the coordinates of the mesh corners. + * @param[out] gradN Array to contain the shape function derivatives for all + * support points at the coordinates of the quadrature point @p q. + * @return The determinant of the parent/physical transformation matrix. + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + static real64 calcGradNWithCorners( real64 const (&coords)[3], + real64 const (&X)[8][3], + real64 ( &gradN )[numNodes][3] ); + + /** + * @brief Calculate the shape functions derivatives wrt the physical + * coordinates. + * @param q Index of the quadrature point. + * @param X Array containing the coordinates of the mesh corners. + * @param stack Variables allocated on the stack as filled by @ref setupStack. + * @param gradN Array to contain the shape function derivatives for all + * support points at the coordinates of the quadrature point @p q. + * @return The determinant of the parent/physical transformation matrix. + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + static real64 calcGradNWithCorners( localIndex const q, + real64 const (&X)[8][3], + StackVariables const & stack, + real64 ( &gradN )[numNodes][3] ); /** * @brief Calculate the integration weights for a quadrature point. @@ -923,6 +968,57 @@ calcGradN( localIndex const q, return calcGradN( q, X, gradN ); } +template< typename GL_BASIS > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradNWithCorners( localIndex const q, + real64 const (&X)[8][3], + real64 (& gradN)[numNodes][3] ) +{ + int qa, qb, qc; + GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc ); + + real64 J[3][3] = {{0}}; + + jacobianTransformation( qa, qb, qc, X, J ); + + real64 const detJ = LvArray::tensorOps::invert< 3 >( J ); + + applyTransformationToParentGradients( q, J, gradN ); + + return detJ; +} +//************************************************************************************************* +template< typename GL_BASIS > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradNWithCorners( real64 const (&coords)[3], + real64 const (&X)[8][3], + real64 (& gradN)[numNodes][3] ) +{ + real64 J[3][3] = {{0}}; + + jacobianTransformation( coords, X, J ); + + real64 const detJ = LvArray::tensorOps::invert< 3 >( J ); + + applyTransformationToParentGradients( coords, J, gradN ); + + return detJ; +} +template< typename GL_BASIS > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >:: +calcGradNWithCorners( localIndex const q, + real64 const (&X)[8][3], + StackVariables const & GEOS_UNUSED_PARAM( stack ), + real64 ( & gradN )[numNodes][3] ) +{ + return calcGradN( q, X, gradN ); +} //************************************************************************************************* #if __GNUC__ #pragma GCC diagnostic push diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp index c6fbf016edd..34539594757 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp @@ -460,13 +460,13 @@ struct PreComputeSourcesAndReceivers sourceCoordinates[isrc][1], sourceCoordinates[isrc][2] }; - real64 xLocal[numNodesPerElem][3]; + real64 xLocal[8][3]; - for( localIndex a=0; a< numNodesPerElem; ++a ) + for( localIndex a = 0; a < 8; ++a ) { - for( localIndex i=0; i<3; ++i ) + for( localIndex i = 0; i < 3; ++i ) { - xLocal[a][i] = baseNodeCoords( elemsToNodes( k, a ), i ); + xLocal[a][i] = baseNodeCoords( elemsToNodes( k, FE_TYPE::meshIndexToLinearIndex3D( a ) ), i ); } } @@ -496,7 +496,7 @@ struct PreComputeSourcesAndReceivers real64 N[numNodesPerElem]; real64 gradN[numNodesPerElem][3]; FE_TYPE::calcN( coordsOnRefElem, N ); - FE_TYPE::calcGradN( coordsOnRefElem, xLocal, gradN ); + FE_TYPE::calcGradNWithCorners( coordsOnRefElem, xLocal, gradN ); R2SymTensor moment = sourceMoment; for( localIndex q=0; q< numNodesPerElem; ++q ) { @@ -593,13 +593,13 @@ struct PreComputeSourcesAndReceivers if( sampleFound && elemGhostRank[k] < 0 ) { real64 coordsOnRefElem[3]{}; - real64 xLocal[numNodesPerElem][3]; + real64 xLocal[8][3]; - for( localIndex a=0; a< numNodesPerElem; ++a ) + for( localIndex a = 0; a < 8; ++a ) { - for( localIndex i=0; i<3; ++i ) + for( localIndex i=0; i < 3; ++i ) { - xLocal[a][i] = baseNodeCoords( elemsToNodes( k, a ), i ); + xLocal[a][i] = baseNodeCoords( elemsToNodes( k, FE_TYPE::meshIndexToLinearIndex3D( a ) ), i ); } } @@ -610,7 +610,7 @@ struct PreComputeSourcesAndReceivers real64 N[numNodesPerElem]; real64 gradN[numNodesPerElem][3]; FE_TYPE::calcN( coordsOnRefElem, N ); - FE_TYPE::calcGradN( coordsOnRefElem, xLocal, gradN ); + FE_TYPE::calcGradNWithCorners( coordsOnRefElem, xLocal, gradN ); for( localIndex a = 0; a < numNodesPerElem; ++a ) { receiverNodeIds[ircv][iSample * numNodesPerElem + a] = elemsToNodes( k, From 4f058cba5a93637aff41e9e9d2a6789532382d65 Mon Sep 17 00:00:00 2001 From: j0405284 Date: Wed, 21 Aug 2024 21:33:34 +0200 Subject: [PATCH 2/6] added corresponding compute jacobian with corners --- .../Qk_Hexahedron_Lagrange_GaussLobatto.hpp | 49 ++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp index 8058ec36baa..74c7a9bffba 100644 --- a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp @@ -539,6 +539,22 @@ class Qk_Hexahedron_Lagrange_GaussLobatto final : public FiniteElementBase static void jacobianTransformation( real64 const (&coords)[3], real64 const (&X)[numNodes][3], real64 ( &J )[3][3] ); + + /** + * @brief Calculates the isoparametric "Jacobian" transformation + * matrix/mapping from the parent space to the physical space at a single point. + * Assumes that the coordinate of high-order nodes are given by trilinear + * interpolation of the mesh corners. + * @param coords The parent coordinates at which to evaluate the shape function value + * @param X Array containing the coordinates of the mesh corners. + * @param J Array to store the Jacobian transformation. + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + static void jacobianTransformationWithCorners( real64 const (&coords)[3], + real64 const (&X)[8][3], + real64 ( &J )[3][3] ); + /** * @brief performs a trilinear interpolation to determine the real-world coordinates of a * vertex @@ -1000,7 +1016,7 @@ Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradNWithCorners( real64 co { real64 J[3][3] = {{0}}; - jacobianTransformation( coords, X, J ); + jacobianTransformationWithCorners( coords, X, J ); real64 const detJ = LvArray::tensorOps::invert< 3 >( J ); @@ -1079,6 +1095,37 @@ jacobianTransformation( real64 const (&coords)[3], }, X, J ); } +template< typename GL_BASIS > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >:: +jacobianTransformationWithCorners( real64 const (&coords)[3], + real64 const (&X)[8][3], + real64 ( & J )[3][3] ) +{ + supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3], + int const nodeIndex, + real64 const (&X)[8][3], + real64 (& J)[3][3] ) + { + int qa, qb, qc; + GL_BASIS::TensorProduct3D::multiIndex( nodeIndex, qa, qb, qc ); + real64 Xnode[3]; + real64 alpha = ( GL_BASIS::parentSupportCoord( qa ) + 1.0 ) / 2.0; + real64 beta = ( GL_BASIS::parentSupportCoord( qb ) + 1.0 ) / 2.0; + real64 gamma = ( GL_BASIS::parentSupportCoord( qc ) + 1.0 ) / 2.0; + trilinearInterp( alpha, beta, gamma, X, Xnode); + for( int i = 0; i < 3; ++i ) + { + for( int j = 0; j < 3; ++j ) + { + J[i][j] = J[i][j] + dNdXi[ j ] * Xnode[i]; + } + } + }, X, J ); +} + template< typename GL_BASIS > GEOS_HOST_DEVICE GEOS_FORCE_INLINE From 2442fd4cc4cd95288b0e6a0f827e853fb4ebfd71 Mon Sep 17 00:00:00 2001 From: j0405284 Date: Wed, 21 Aug 2024 21:37:59 +0200 Subject: [PATCH 3/6] uncrustify --- .../elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp index 74c7a9bffba..bb49c000c71 100644 --- a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp @@ -1115,7 +1115,7 @@ jacobianTransformationWithCorners( real64 const (&coords)[3], real64 alpha = ( GL_BASIS::parentSupportCoord( qa ) + 1.0 ) / 2.0; real64 beta = ( GL_BASIS::parentSupportCoord( qb ) + 1.0 ) / 2.0; real64 gamma = ( GL_BASIS::parentSupportCoord( qc ) + 1.0 ) / 2.0; - trilinearInterp( alpha, beta, gamma, X, Xnode); + trilinearInterp( alpha, beta, gamma, X, Xnode ); for( int i = 0; i < 3; ++i ) { for( int j = 0; j < 3; ++j ) From 064e9e5610e9d928f1d49b7c43c9e0ec7f70ed64 Mon Sep 17 00:00:00 2001 From: j0405284 Date: Thu, 22 Aug 2024 18:07:12 +0200 Subject: [PATCH 4/6] added correction for ghost nodes --- .../shared/PrecomputeSourcesAndReceiversKernel.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp index 34539594757..741a2a80c78 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp @@ -466,7 +466,7 @@ struct PreComputeSourcesAndReceivers { for( localIndex i = 0; i < 3; ++i ) { - xLocal[a][i] = baseNodeCoords( elemsToNodes( k, FE_TYPE::meshIndexToLinearIndex3D( a ) ), i ); + xLocal[a][i] = baseNodeCoords( baseElemsToNodes( k, a ), i ); } } @@ -599,7 +599,7 @@ struct PreComputeSourcesAndReceivers { for( localIndex i=0; i < 3; ++i ) { - xLocal[a][i] = baseNodeCoords( elemsToNodes( k, FE_TYPE::meshIndexToLinearIndex3D( a ) ), i ); + xLocal[a][i] = baseNodeCoords( baseElemsToNodes( k, a ), i ); } } From d3d9f5ef8109b498fd6750168ed141fbcdba3c9f Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Thu, 29 Aug 2024 09:20:28 -0700 Subject: [PATCH 5/6] Re-enable wave tests --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index a00377b2bbc..3fa9cbb2698 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -4,4 +4,4 @@ baselines: allow_fail: all: '' - streak: pennyShapedToughnessDominated_smoke_01,pennyShapedViscosityDominated_smoke_01,pknViscosityDominated_smoke_01,SeismicityRate_poromechanics_1d_smoke_01,SeismicityRate_poromechanics_1d_smoke_02,elas3D_Q3_abc_smoke_01,elas3D_Q3_abc_smoke_08,elas3D_Q3_abc_fs_smoke_01,elas3D_Q3_abc_fs_smoke_08 + streak: pennyShapedToughnessDominated_smoke_01,pennyShapedViscosityDominated_smoke_01,pknViscosityDominated_smoke_01,SeismicityRate_poromechanics_1d_smoke_01,SeismicityRate_poromechanics_1d_smoke_02 From 516975eee22d2cbc7b8e190a761864ea73e1ad68 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Sun, 1 Sep 2024 16:03:23 -0700 Subject: [PATCH 6/6] Update BASELINE_NOTES.md --- BASELINE_NOTES.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index d1d315d5658..13eca60093d 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ 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 #3294 (2024-09-01) +====================== +Re-enable enforcement of wave propagation integrated test pass. + PR #3300 (2024-08-28) ====================== Re-enable floating point exceptions. Rebaseline due to minor changing default value of maxRelativeCompDensChange from 1.7976931348623157e+308 to 1.7976931348623157e+208.