Skip to content

Commit

Permalink
fix: Bug with elastic wave propagator at higher orders (#3294)
Browse files Browse the repository at this point in the history 
* fixed a bug due to PR 3202 concerning the elastic solver at higher orders

* added corresponding compute jacobian with corners

* added correction for ghost nodes

* Re-enable wave tests

---------

Co-authored-by: acitrain <[email protected]>
  • Loading branch information
@sframba @acitrain
sframba and acitrain authored Sep 2, 2024
1 parent 5f6d710 commit a38d1ee
Show file tree
Hide file tree
Showing 4 changed files with 158 additions and 11 deletions.
2 changes: 1 addition & 1 deletion .integrated_tests.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
4 changes: 4 additions & 0 deletions BASELINE_NOTES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
143 changes: 143 additions & 0 deletions src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -494,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
Expand Down Expand Up @@ -923,6 +984,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}};

jacobianTransformationWithCorners( 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
Expand Down Expand Up @@ -983,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
Expand Down
20 changes: 10 additions & 10 deletions ...eComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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( baseElemsToNodes( k, a ), i );
}
}

Expand Down Expand Up @@ -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 )
{
Expand Down Expand Up @@ -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( baseElemsToNodes( k, a ), i );
}
}

Expand All @@ -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,
Expand Down

0 comments on commit a38d1ee

Please sign in to comment.