AcousticWaveEquationDGKernel.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2020 TotalEnergies
 * Copyright (c) 2019-     GEOS Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICWAVEEQUATIONDGKERNEL_HPP_
#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICWAVEEQUATIONDGKERNEL_HPP_
 
#include "finiteElement/kernelInterface/KernelBase.hpp"
#include "denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp"
 
namespace geos
{
 
namespace AcousticWaveEquationDGKernels
{
 
struct PrecomputeSourceAndReceiverKernel
{
  using EXEC_POLICY = parallelDevicePolicy< >;
 
  template< typename EXEC_POLICY, typename FE_TYPE >
  static void
  launch( localIndex const size,
          localIndex const regionIndex,
          ArrayOfArraysView< localIndex const > const baseFacesToNodes,
          arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const baseNodeCoords,
          arrayView1d< globalIndex const > const baseNodeLocalToGlobal,
          arrayView1d< globalIndex const > const elementLocalToGlobal,
          ArrayOfArraysView< localIndex const > const baseNodesToElements,
          arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes,
          arrayView1d< integer const > const elemGhostRank,
          arrayView2d< localIndex const > const elemsToFaces,
          arrayView2d< real64 const > const & elemCenter,
          arrayView2d< real64 const > const sourceCoordinates,
          arrayView1d< localIndex > const sourceIsAccessible,
          arrayView1d< localIndex > const sourceElem,
          arrayView2d< real64 > const sourceConstants,
          arrayView1d< localIndex > const sourceRegion,
          arrayView2d< real64 const > const receiverCoordinates,
          arrayView1d< localIndex > const receiverIsLocal,
          arrayView1d< localIndex > const receiverElem,
          arrayView2d< real64 > const receiverConstants,
          arrayView1d< localIndex > const receiverRegion )
  {
 
    forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
    {
      real64 const center[3] = { elemCenter[k][0],
                                 elemCenter[k][1],
                                 elemCenter[k][2] };
 
      // Step 1: locate the sources, and precompute the source term
 
      for( localIndex isrc = 0; isrc < sourceCoordinates.size( 0 ); ++isrc )
      {
        if( sourceIsAccessible[isrc] == 0 )
        {
          real64 const coords[3] = { sourceCoordinates[isrc][0],
                                     sourceCoordinates[isrc][1],
                                     sourceCoordinates[isrc][2] };
          bool const sourceFound =
            computationalGeometry::isPointInsideConvexPolyhedronRobust( k,
                                                                        baseNodeCoords,
                                                                        elemsToFaces,
                                                                        baseFacesToNodes,
                                                                        baseNodesToElements,
                                                                        baseNodeLocalToGlobal,
                                                                        elementLocalToGlobal,
                                                                        center,
                                                                        coords );
 
          if( sourceFound )
          {
            sourceIsAccessible[isrc] = 1;
            sourceElem[isrc] = k;
            sourceRegion[isrc] = regionIndex;
            real64 Ntest[FE_TYPE::numNodes];
 
 
 
            constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
            real64 xLocal[4][3];
            for( localIndex a=0; a<4; ++a )
            {
              for( localIndex i=0; i<3; ++i )
              {
                xLocal[a][i] = baseNodeCoords( baseElemsToNodes( k, a ), i );
              }
            }
 
            FE_TYPE::calcN( xLocal, coords, Ntest );
 
            for( localIndex a = 0; a < numNodesPerElem; ++a )
            {
              sourceConstants[isrc][a] = Ntest[a];
            }
          }
        }
      } // end loop over all sources
 
      // Step 2: locate the receivers, and precompute the receiver term
 
      for( localIndex ircv = 0; ircv < receiverCoordinates.size( 0 ); ++ircv )
      {
        if( receiverIsLocal[ircv] == 0 )
        {
          real64 const coords[3] = { receiverCoordinates[ircv][0],
                                     receiverCoordinates[ircv][1],
                                     receiverCoordinates[ircv][2] };
 
          bool const receiverFound =
            computationalGeometry::isPointInsideConvexPolyhedronRobust( k,
                                                                        baseNodeCoords,
                                                                        elemsToFaces,
                                                                        baseFacesToNodes,
                                                                        baseNodesToElements,
                                                                        baseNodeLocalToGlobal,
                                                                        elementLocalToGlobal,
                                                                        center,
                                                                        coords );
 
          if( receiverFound && elemGhostRank[k] < 0 )
          {
            receiverIsLocal[ircv] = 1;
            receiverElem[ircv] = k;
            receiverRegion[ircv] = regionIndex;
 
 
            real64 Ntest[FE_TYPE::numNodes];
            constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
 
            real64 xLocal[4][3];
            for( localIndex a=0; a< 4; ++a )
            {
              for( localIndex i=0; i<3; ++i )
              {
                xLocal[a][i] = baseNodeCoords( baseElemsToNodes( k, a ), i );
              }
            }
            FE_TYPE::calcN( xLocal, coords, Ntest );
 
            for( localIndex a = 0; a < numNodesPerElem; ++a )
            {
              receiverConstants[ircv][a] = Ntest[a];
            }
          }
        }
      } // end loop over receivers
 
    } );
 
  }
};
 
struct PrecomputeNeighborhoodKernel
{
 
  using EXEC_POLICY = parallelDevicePolicy< >;
 
  template< typename EXEC_POLICY, typename FE_TYPE >
  static void
  launch( localIndex const size,
          arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes,
          arrayView2d< localIndex const > const & elemsToFaces,
          arrayView2d< localIndex const > const & facesToElems,
          ArrayOfArraysView< localIndex const > const & facesToNodes,
          arrayView1d< localIndex const > const & freeSurfaceFaceIndicator,
          arrayView2d< localIndex > const & elemsToOpposite,
          arrayView2d< integer > const & elemsToOppositePermutation )
  {
    forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k1 )
    {
      localIndex vertices[ 4 ] = { elemsToNodes( k1, 0 ), elemsToNodes( k1, 1 ), elemsToNodes( k1, 2 ), elemsToNodes( k1, 3 ) };
      for( int i = 0; i < 4; i++ )
      {
        localIndex k1OrderedVertices[ 3 ];
        localIndex f = elemsToFaces( k1, i );
        localIndex faceVertices[ 3 ] = { facesToNodes( f, 0 ), facesToNodes( f, 1 ), facesToNodes( f, 2 ) };
        // find neighboring element, if any
        localIndex k2 = facesToElems( f, 0 );
        if( k2 == k1 )
        {
          k2 = facesToElems( f, 1 );
        }
        // find opposite vertex in first element
        int o1 = -1;
        int indexo1= -1;
        int vertex = -1;
        int count = 0;
        for( localIndex k=0; k< 4; ++k )
        {
          vertex = vertices[k];
          bool found = false;
          for( int j = 0; j < 3; j++ )
          {
            if( vertex == faceVertices[ j ] )
            {
              found = true;
              break;
            }
          }
          if( !found )
          {
            o1 = vertex;
            indexo1=k;
          }
          else
          {
            k1OrderedVertices[ count++ ] = vertex;
          }
        }
 
        GEOS_ERROR_IF( o1 < 0, "Topological error in mesh: a face and its adjacent element share all vertices." );
        if( k2 < 0 )
        {
          // boundary element, either free surface, or absorbing boundary
          elemsToOpposite( k1, indexo1 ) = freeSurfaceFaceIndicator( f ) == 1 ? -2 : -1;
          elemsToOppositePermutation( k1, indexo1 ) = 0;
        }
        else
        {
          elemsToOpposite( k1, indexo1 ) = k2;
          localIndex oppositeElemVertices[ 4 ] = { elemsToNodes( k2, 0 ), elemsToNodes( k2, 1 ), elemsToNodes( k2, 2 ), elemsToNodes( k2, 3 ) };
          // compute permutation
          integer permutation = 0;
          int c = 1;
          for( localIndex k2Vertex : oppositeElemVertices )
          {
            int position = -1;
            for( int j = 0; j < 3; j++ )
            {
              if( k1OrderedVertices[ j ] == k2Vertex )
              {
                position = j;
                break;
              }
            }
            permutation = permutation + c * ( position + 1 );
            c = c * 4;
          }
          elemsToOppositePermutation( k1, indexo1 ) = permutation;
        }
      }
    } );
  }
};
 
struct PrecomputePenaltyGeomKernel
{
  using EXEC_POLICY = parallelDevicePolicy< >;
 
  template< typename EXEC_POLICY >
  static void
  launch( localIndex const size,
          arrayView2d< real32 const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
          arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes,
          arrayView1d< real32 > const & characteristicSize )
  {
    forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
    {
      characteristicSize[ k ] = WaveSolverUtils::computeReferenceLengthForPenalty( elemsToNodes[ k ], nodeCoords );
    } );
  }
};
 
 
//TODO : in the future if needed for boundary condition
//struct PrecomputeMassDampingKernel
//{
//  using EXEC_POLICY = parallelDevicePolicy< >;
//
//  /**
//   * @brief Launches the precomputation of the inverse of the reference mass matrix in the bulk, as well as
//   *   the inverse mass + damping term for the boundary elements.
//   * @tparam EXEC_POLICY execution policy
//   * @tparam FE_TYPE finite element type
//   * @param[in] size the number of elements
//   * @param[in] nodeCoords coordinates of the nodes
//   * @param[in] elemsToNodes map from element to nodes
//   * @param[in] elemToOpposite DG element-to-opposite map
//   * @param[out] referenceInvMassMatrix computed M^{-1} for the reference element
//   * @param[out] boundaryInvMassPlusDamping (M + dt/2 D)^{-1} for boundary elements
//   * @param[in] dt time-step
//   */
//  template< typename EXEC_POLICY, typename ATOMIC_POLICY, typename FE_TYPE >
//  static void
//  launch( localIndex const size,
//          arrayView2d< real32 const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
//          arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes,
//          arrayView2d< localIndex const > const & elemsToOpposite,
//          array2d< real64 > & referenceInvMassMatrix,
//          array3d< real64 > & boundaryInvMassPlusDamping,
//          real64 const dt )
//  {
//    // Precompute reference mass matrix for non-boundary elements
//    referenceInvMassMatrix.resizeDimension< 0, 1 >( FE_TYPE::numNodes, FE_TYPE::numNodes );
//    array2d< real64 > massMatrix;
//    massMatrix.resize( FE_TYPE::numNodes, FE_TYPE::numNodes );
//    massMatrix.zero();
//    FE_TYPE::computeReferenceMassMatrix( massMatrix );
//    BlasLapackLA::matrixInverse( massMatrix, referenceInvMassMatrix );
//    // Precompute local mass + damping matrix on the boundary elements
//    localIndex nAbsBdryElems = 0;
//    forAll< EXEC_POLICY >( size, [&] ( localIndex const k )
//    {
//      bool bdry = false;
//      for( int i = 0; i < 4; i++ )
//      {
//        if( elemsToOpposite[ k ][ i ] == -1 )
//        {
//          bdry = true;
//          break;
//        }
//      }
//      if( bdry )
//      {
//        RAJA::atomicInc< ATOMIC_POLICY >( &nAbsBdryElems );
//      }
//    } );
//
//    boundaryInvMassPlusDamping.resizeDimension< 0, 1, 2 >( nAbsBdryElems, FE_TYPE::numNodes, FE_TYPE::numNodes );
//    forAll< EXEC_POLICY >( size, [&] ( localIndex const k )
//    {} );
//  }
//};
 
 
 
struct PressureComputationKernel
{
  using EXEC_POLICY = parallelDevicePolicy< >;
 
  template< typename FE_TYPE, typename EXEC_POLICY, typename ATOMIC_POLICY >
  static void
  pressureComputation( localIndex const size,
                       localIndex const regionIndex,
                       arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
                       arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes,
                       arrayView1d< real32 const > const velocity,
                       arrayView1d< real32 const > const density,
                       arrayView2d< real32 > const p_n,
                       arrayView2d< real32 const > const p_nm1,
                       arrayView2d< localIndex > const & elemsToOpposite,
                       arrayView2d< integer > const & elemsToOppositePermutation,
                       arrayView2d< real64 >referenceInvMassMatrix,
                       arrayView1d< real32 const > const characteristicSize,
                       arrayView2d< real64 const > const sourceConstants,
                       arrayView1d< localIndex const > const sourceIsAccessible,
                       arrayView1d< localIndex const > const sourceElem,
                       arrayView1d< localIndex const > const sourceRegion,
                       real64 const dt,
                       real64 const time_n,
                       real32 const timeSourceFrequency,
                       real32 const timeSourceDelay,
                       localIndex const rickerOrder,
                       bool const useSourceWaveletTables,
                       arrayView1d< TableFunction::KernelWrapper const > const sourceWaveletTableWrappers,
                       arrayView2d< real32 > const p_np1 )
 
  {
 
    real64 const rickerValue = useSourceWaveletTables ? 0 : WaveSolverUtils::evaluateRicker( time_n, timeSourceFrequency, timeSourceDelay, rickerOrder );
 
    forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
    {
 
      int const order = FE_TYPE::order;
 
      // Coefficient for penalisation
      int const alpha = 12.0*((order*(order+1)*(order+2))/6);
 
      real64 const dt2 = pow( dt, 2 );
 
      constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
 
      real64 pTemp[numNodesPerElem] = {0.0};
      real64 flowx[numNodesPerElem] = {0.0};
 
      real64 xLocal[4][3];
      for( localIndex a=0; a< 4; ++a )
      {
        for( localIndex i=0; i<3; ++i )
        {
          xLocal[a][i] = X( elemsToNodes( k, a ), i );
        }
      }
 
      //Compute 1/(c2*density) term and compute the determinant of the jacobian
      real64 const C2 = pow( velocity[k], 2 );
 
      real64 const invC2 = 1.0/(C2*density[k]);
 
      real64 const det = LvArray::math::abs( FE_TYPE::jacobianDeterminant( xLocal ));
 
      //Multiply p_{n } by 2*invC2*Mass and p_{n-1} by -invC2*Mass
      FE_TYPE::computeMassTerm( xLocal, [&] ( const int i, const int j, const real64 val )
      {
        pTemp[i] += 2.0*invC2*val*p_n[k][j];
        pTemp[i] -= val*invC2*p_nm1[k][j];
      } );
 
      //First stiffness part (volume)
      FE_TYPE::computeStiffnessTerm( xLocal, [&] ( const int i, const int j, real64 val )
      {
        flowx[i] -= val*p_n[k][j];
      } );
 
 
      //Second stiffness part (surface)
      FE_TYPE::computeSurfaceTerms( xLocal, [&] ( const int c1, const int c2, const int f1, const int, const int, const int, const int i2, const int j2, const int k2, real64 const val )
      {
 
        // We take the neighbour element
        const localIndex elemNeigh = elemsToOpposite( k, f1 );
 
 
        if( elemNeigh >= 0 )
        {
 
 
          // We compute the permutation of the opposite element, used to find the equivalent of c2 on the neighbour element
          const int perm = elemsToOppositePermutation( k, f1 );
 
 
          const int p1 = perm%4-1;
          const int p2 = (perm/4)%4-1;
          const int p3 = (perm/16)%4-1;
 
          const int l2 = order-i2-j2-k2;
          const int Indices[3] = {i2, j2, k2};
 
          const int ii2 = p1 < 0 ? l2 : Indices[p1];
          const int jj2 = p2 < 0 ? l2 : Indices[p2];
          const int kk2 = p3 < 0 ? l2 : Indices[p3];
 
          const int neighDof2 = FE_TYPE::dofIndex( ii2, jj2, kk2 );
 
          //Now that we have the neighbor dof, we compute the penalisation term. The extra tyerm 12 is used as a safety for the value
          // of coefficient alpha= 1/hmin
 
 
          real64 const val1 = (alpha/((LvArray::math::min( characteristicSize[k], characteristicSize[elemNeigh] ))))*val*p_n[k][c2];
          real64 const val2 = (alpha/((LvArray::math::min( characteristicSize[k], characteristicSize[elemNeigh] ))))*val*p_n[elemNeigh][neighDof2];
 
          flowx[c1] += -val1+val2;
 
        }
 
      },
                                    [&] ( const int c1, const int c2, const int f1, const int fNeigh, const int i1, const int j1, const int k1, const int i2, const int j2, const int k2,
                                          real64 const val )
      {
 
 
        // This part is for the flux term
        //We take the neighbour element
        const int elemNeigh = elemsToOpposite( k, f1 );
 
        if( elemNeigh >= 0 )
        {
 
          const int perm = elemsToOppositePermutation( k, f1 );
 
          const int p1 = perm%4-1;
          const int p2 = (perm/4)%4-1;
          const int p3 = (perm/16)%4-1;
          const int p4 = (perm/64)%4-1;
 
          // Again we take the permutation of the opposite element, used to find the equivalent of c2 on the neighbour element and c1 this
          // time
 
          const int Indices[3] = {i2, j2, k2};
 
          const int l2 = order-i2-j2-k2;
 
          const int ii2 = p1 < 0 ? l2 : Indices[p1];
          const int jj2 = p2 < 0 ? l2 : Indices[p2];
          const int kk2 = p3 < 0 ? l2 : Indices[p3];
 
          const int neighDof2 = FE_TYPE::dofIndex( ii2, jj2, kk2 );
 
          const int IndicesTranspose[3] = {i1, j1, k1};
 
          const int l1 = order-i1-j1-k1;
 
          const int ii1 = p1 < 0 ? l1 : IndicesTranspose[p1];
          const int jj1 = p2 < 0 ? l1 : IndicesTranspose[p2];
          const int kk1 = p3 < 0 ? l1 : IndicesTranspose[p3];
 
          const int neighDof = FE_TYPE::dofIndex( ii1, jj1, kk1 );
 
          // Here the goal is to compute the factor of correction for the flux term
          // We need to get the local coordinates of the element and the neighbour element on a specific order
 
          //We first assign to elemneigh the three value from xLocal which are common to the two elements
          real64 xLocalNeigh[4][3];
 
          int x1 = (f1+1)%4;
          int x2 = (f1+2)%4;
          int x3 = (f1+3)%4;
 
          if( x2 == fNeigh )
          {
            int temp = x2;
            x2 = x3;
            x3 = temp;
          }
 
          if( x1 == fNeigh )
          {
            int temp = x1;
            x1 = x3;
            x3 = temp;
          }
 
          int o1 = fNeigh < 0 ? f1 : fNeigh;
          int o2 = f1;
 
          xLocalNeigh[x1][0] = xLocal[x1][0];
          xLocalNeigh[x1][1] = xLocal[x1][1];
          xLocalNeigh[x1][2] = xLocal[x1][2];
          xLocalNeigh[x2][0] = xLocal[x2][0];
          xLocalNeigh[x2][1] = xLocal[x2][1];
          xLocalNeigh[x2][2] = xLocal[x2][2];
          xLocalNeigh[x3][0] = xLocal[x3][0];
          xLocalNeigh[x3][1] = xLocal[x3][1];
          xLocalNeigh[x3][2] = xLocal[x3][2];
 
          // Using the permutation we can find the last point of the neighbour element
 
          localIndex IndexPerm[4] = {p1, p2, p3, p4};
 
          for( localIndex i = 0; i <4; ++i )
          {
            if( IndexPerm[i] ==-1 )
            {
              xLocalNeigh[f1][0] = X[ elemsToNodes( elemNeigh, i )][0];
              xLocalNeigh[f1][1] = X[ elemsToNodes( elemNeigh, i )][1];
              xLocalNeigh[f1][2] = X[ elemsToNodes( elemNeigh, i )][2];
            }
 
          }
 
          //We compute the correction factor and the value of the flux term
          real64 corrLocal = FE_TYPE::computeFluxDerivativeFactor( xLocal, x1, x2, o1, o2 );
          real64 corrNeigh = FE_TYPE::computeFluxDerivativeFactor( xLocalNeigh, x1, x2, o1, o2 );
 
          real64 const val1 = 0.5*val*p_n[k][c2]*corrLocal;
          real64 const val2 = 0.5*val*p_n[elemNeigh][neighDof2]*corrLocal;
//
          real64 const val3 = 0.5*val*p_n[k][c1]*corrLocal;
          real64 const val4 = 0.5*val*p_n[elemNeigh][neighDof]*corrNeigh;
 
          flowx[c1] += val1-val2;
          flowx[c2] += val3-val4;
 
        }
      } );
 
      //Add time dependency
      for( localIndex i = 0; i < numNodesPerElem; i++ )
      {
        pTemp[i] += dt2*flowx[i];
      }
 
      //Source Injection
      for( localIndex isrc = 0; isrc < sourceConstants.size( 0 ); ++isrc )
      {
        if( sourceIsAccessible[isrc] == 1 )
        {
          if( sourceElem[isrc]==k && sourceRegion[isrc] == regionIndex )
          {
 
            real64 const srcValue = useSourceWaveletTables ? sourceWaveletTableWrappers[ isrc ].compute( &time_n ) : rickerValue;
            for( localIndex i = 0; i < numNodesPerElem; ++i )
            {
              pTemp[i]+= dt2*(sourceConstants[isrc][i]*srcValue);
            }
          }
        }
      }
 
      //Multiplication by the inverse mass matrix
      for( localIndex i = 0; i < numNodesPerElem; ++i )
      {
        real64 fx=0.0;
        for( localIndex j = 0; j < numNodesPerElem; ++j )
        {
          fx+= referenceInvMassMatrix( i, j )*pTemp[j];
        }
        p_np1[k][i]=(C2*density[k]*fx)/det;
      }
 
    } );
 
 
 
  }
 
  // FE_TYPE const & m_finiteElement;
 
 
};
 
} // namespace AcousticWaveEquationDGKernels
 
} // namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICWAVEEQUATIONDGKERNEL_HPP_