WaveSolverUtils.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERUTILS_HPP_
#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERUTILS_HPP_
 
#include "mesh/utilities/ComputationalGeometry.hpp"
#include "fileIO/Outputs/OutputBase.hpp"
#include "LvArray/src/tensorOps.hpp"
 
namespace geos
{
 
struct WaveSolverUtils
{
  static constexpr real64 epsilonLoc = 1e-8;
 
  using EXEC_POLICY = parallelDevicePolicy<>;
  using wsCoordType = real32;
 
  enum class DASType : integer
  {
    none,                
    dipole,              
    strainIntegration,   
  };
 
  enum class AttenuationType : integer
  {
    none,   
    sls,    
  };
 
  GEOS_HOST_DEVICE
  static real32 evaluateRicker( real64 const time_n, real32 const f0, real32 const t0, localIndex const order )
  {
    real32 const delay = t0 > 0 ? t0 : 1 / f0;
    real32 pulse = 0.0;
    real32 const alpha = -pow( f0 * M_PI, 2 );
    real32 const time_d = time_n - delay;
    real32 const gaussian = exp( alpha * pow( time_d, 2 ));
    localIndex const sgn = pow( -1, order + 1 );
 
    switch( order )
    {
      case 0:
        pulse = sgn * gaussian;
        break;
      case 1:
        pulse = sgn * (2 * alpha * time_d) * gaussian;
        break;
      case 2:
        pulse = sgn * (2 * alpha + 4 * pow( alpha, 2 ) * pow( time_d, 2 )) * gaussian;
        break;
      case 3:
        pulse = sgn * (12 * pow( alpha, 2 ) * time_d + 8 * pow( alpha, 3 ) * pow( time_d, 3 )) * gaussian;
        break;
      case 4:
        pulse = sgn * (12 * pow( alpha, 2 ) + 48 * pow( alpha, 3 ) * pow( time_d, 2 ) + 16 * pow( alpha, 4 ) * pow( time_d, 4 )) * gaussian;
        break;
      default:
        GEOS_ERROR( "This option is not supported yet, rickerOrder must be in range {0:4}" );
    }
 
    return pulse;
  }
 
  static void initTrace( char const *prefix,
                         string const & name,
                         bool const outputSeismoTrace,
                         localIndex const nReceivers,
                         arrayView1d< localIndex const > const receiverIsLocal )
  {
    if( !outputSeismoTrace )
      return;
 
    string const outputDir = OutputBase::getOutputDirectory();
    RAJA::ReduceSum< ReducePolicy< serialPolicy >, localIndex > count( 0 );
 
    forAll< serialPolicy >( nReceivers, [=]( localIndex const ircv )
    {
      if( receiverIsLocal[ircv] == 1 )
      {
        count += 1;
        string const fn = joinPath( outputDir, GEOS_FMT( "{}_{}_{:03}.txt", prefix, name, ircv ) );
        std::ofstream f( fn, std::ios::out | std::ios::trunc );
      }
    } );
 
    localIndex const total = MpiWrapper::sum( count.get());
    GEOS_ERROR_IF( nReceivers != total, GEOS_FMT( ": Invalid distribution of receivers: nReceivers={} != MPI::sum={}.", nReceivers, total ));
  }
 
  static void writeSeismoTraceVector( char const *prefix,
                                      string const & name,
                                      bool const outputSeismoTrace,
                                      localIndex const nReceivers,
                                      arrayView1d< localIndex const > const receiverIsLocal,
                                      localIndex const nsamplesSeismoTrace,
                                      arrayView2d< real32 const > const varAtReceiversx,
                                      arrayView2d< real32 const > const varAtReceiversy,
                                      arrayView2d< real32 const > const varAtReceiversz )
  {
    writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversx );
    writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversy );
    writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversz );
  }
 
  static void writeSeismoTrace( char const *prefix,
                                string const & name,
                                bool const outputSeismoTrace,
                                localIndex const nReceivers,
                                arrayView1d< localIndex const > const receiverIsLocal,
                                localIndex const nsamplesSeismoTrace,
                                arrayView2d< real32 const > const varAtReceivers )
  {
    if( !outputSeismoTrace )
      return;
 
    string const outputDir = OutputBase::getOutputDirectory();
    forAll< serialPolicy >( nReceivers, [=]( localIndex const ircv )
    {
      if( receiverIsLocal[ircv] == 1 )
      {
        string const fn = joinPath( outputDir, GEOS_FMT( "{}_{}_{:03}.txt", prefix, name, ircv ) );
        std::ofstream f( fn, std::ios::app );
        if( f )
        {
          GEOS_LOG_RANK( GEOS_FMT( "Append to seismo trace file {}", fn ) );
          for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample )
          {
            // index - time - value
            f << iSample << " " << varAtReceivers[iSample][nReceivers] << " " << varAtReceivers[iSample][ircv] << std::endl;
          }
          f.close();
        }
        else
        {
          GEOS_WARNING( GEOS_FMT( "Failed to open output file {}", fn ) );
        }
      }
    } );
  }
 
  static void computeSeismoTrace( real64 const time_n,
                                  real64 const dt,
                                  real64 const timeSeismo,
                                  localIndex const iSeismo,
                                  arrayView2d< localIndex const > const receiverNodeIds,
                                  arrayView2d< real64 const > const receiverConstants,
                                  arrayView1d< localIndex const > const receiverIsLocal,
                                  arrayView1d< real32 const > const var_np1,
                                  arrayView1d< real32 const > const var_n,
                                  arrayView2d< real32 > varAtReceivers,
                                  arrayView1d< real32 > coeffs = {},
                                  bool add = false )
  {
    real64 const time_np1 = time_n + dt;
 
    real32 const a1 = LvArray::math::abs( dt ) < epsilonLoc ? 1.0 : (time_np1 - timeSeismo) / dt;
    real32 const a2 = 1.0 - a1;
 
    localIndex const nReceivers = receiverConstants.size( 0 );
    forAll< EXEC_POLICY >( nReceivers, [=] GEOS_HOST_DEVICE ( localIndex const ircv )
    {
      if( receiverIsLocal[ircv] > 0 )
      {
        real32 vtmp_np1 = 0.0, vtmp_n = 0.0;
        for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode )
        {
          if( receiverNodeIds( ircv, inode ) >= 0 )
          {
            vtmp_np1 += var_np1[receiverNodeIds( ircv, inode )] * receiverConstants( ircv, inode );
            vtmp_n += var_n[receiverNodeIds( ircv, inode )] * receiverConstants( ircv, inode );
          }
        }
        // linear interpolation between the pressure value at time_n and time_{n+1}
        real32 receiverCoeff = coeffs.size( 0 ) == 0 ? 1.0 : coeffs( ircv );
        if( add )
        {
          varAtReceivers( iSeismo, ircv ) += receiverCoeff * ( a1 * vtmp_n + a2 * vtmp_np1 );
        }
        else
        {
          varAtReceivers( iSeismo, ircv ) = receiverCoeff * ( a1 * vtmp_n + a2 * vtmp_np1 );
        }
        // NOTE: varAtReceivers has size(1) = numReceiversGlobal + 1, this does not OOB
        // left in the forAll loop for sync issues since the following does not depend on `ircv`
        varAtReceivers( iSeismo, nReceivers ) = a1 * time_n + a2 * time_np1;
      }
    } );
  }
 
  static void compute2dVariableSeismoTrace( real64 const time_n,
                                            real64 const dt,
                                            localIndex const regionIndex,
                                            arrayView1d< localIndex const > const receiverRegion,
                                            real64 const timeSeismo,
                                            localIndex const iSeismo,
                                            arrayView1d< localIndex const > const receiverElem,
                                            arrayView2d< real64 const > const receiverConstants,
                                            arrayView1d< localIndex const > const receiverIsLocal,
                                            arrayView2d< real32 const > const var_np1,
                                            arrayView2d< real32 const > const var_n,
                                            arrayView2d< real32 > varAtReceivers )
  {
    real64 const time_np1 = time_n + dt;
 
    real32 const a1 = dt < epsilonLoc ? 1.0 : (time_np1 - timeSeismo) / dt;
    real32 const a2 = 1.0 - a1;
 
    localIndex const nReceivers = receiverConstants.size( 0 );
 
    forAll< EXEC_POLICY >( nReceivers, [=] GEOS_HOST_DEVICE ( localIndex const ircv )
    {
      if( receiverIsLocal[ircv] == 1 )
      {
        if( receiverRegion[ircv] == regionIndex )
        {
          real32 vtmp_np1 = 0.0, vtmp_n = 0.0;
          for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode )
          {
            vtmp_np1 += var_np1( receiverElem[ircv], inode ) * receiverConstants( ircv, inode );
            vtmp_n += var_n( receiverElem[ircv], inode ) * receiverConstants( ircv, inode );
          }
          // linear interpolation between the pressure value at time_n and time_{n+1}
          varAtReceivers( iSeismo, ircv ) = a1 * vtmp_n + a2 * vtmp_np1;
          // NOTE: varAtReceivers has size(1) = numReceiversGlobal + 1, this does not OOB
          // left in the forAll loop for sync issues since the following does not depend on `ircv`
          varAtReceivers( iSeismo, nReceivers ) = a1 * time_n + a2 * time_np1;
        }
      }
    } );
  }
 
  GEOS_HOST_DEVICE
  static bool
  locateSourceElement( real64 const numFacesPerElem,
                       real64 const (&elemCenter)[3],
                       arrayView2d< real64 const > const faceNormal,
                       arrayView2d< real64 const > const faceCenter,
                       arraySlice1d< localIndex const > const elemsToFaces,
                       real64 const (&coords)[3] )
  {
    // Loop over the element faces
    real64 tmpVector[3]{};
    for( localIndex kfe = 0; kfe < numFacesPerElem; ++kfe )
    {
 
      localIndex const iface = elemsToFaces[kfe];
      real64 faceCenterOnFace[3] = {faceCenter[iface][0],
                                    faceCenter[iface][1],
                                    faceCenter[iface][2]};
      real64 faceNormalOnFace[3] = {faceNormal[iface][0],
                                    faceNormal[iface][1],
                                    faceNormal[iface][2]};
 
      // Test to make sure if the normal is outwardly directed
      LvArray::tensorOps::copy< 3 >( tmpVector, faceCenterOnFace );
      LvArray::tensorOps::subtract< 3 >( tmpVector, elemCenter );
      if( LvArray::tensorOps::AiBi< 3 >( tmpVector, faceNormalOnFace ) < 0.0 )
      {
        LvArray::tensorOps::scale< 3 >( faceNormalOnFace, -1 );
      }
 
      // compute the vector face center to query point
      LvArray::tensorOps::subtract< 3 >( faceCenterOnFace, coords );
      localIndex const s = computationalGeometry::sign( LvArray::tensorOps::AiBi< 3 >( faceNormalOnFace, faceCenterOnFace ));
 
      // all dot products should be non-negative (we enforce outward normals)
      if( s < 0 )
        return false;
    }
    return true;
  }
 
  template< typename FE_TYPE >
  GEOS_HOST_DEVICE static void
  computeCoordinatesOnReferenceElement( real64 const (&coords)[3],
                                        arraySlice1d< localIndex const, cells::NODE_MAP_USD - 1 > const elemsToNodes,
                                        arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
                                        real64 (& coordsOnRefElem)[3] )
  {
    // only the eight corners of the mesh cell are needed to compute the Jacobian
    real64 xLocal[8][3]{};
    for( localIndex a = 0; a < 8; ++a )
    {
      LvArray::tensorOps::copy< 3 >( xLocal[a], nodeCoords[elemsToNodes[a]] );
    }
    // coordsOnRefElem = invJ*(coords-coordsNode_0)
    real64 invJ[3][3]{};
    FE_TYPE::invJacobianTransformation( 0, xLocal, invJ );
    for( localIndex i = 0; i < 3; ++i )
    {
      // init at (-1,-1,-1) as the origin of the referential elem
      coordsOnRefElem[i] = -1.0;
      for( localIndex j = 0; j < 3; ++j )
      {
        coordsOnRefElem[i] += invJ[i][j] * (coords[j] - xLocal[0][j]);
      }
    }
  }
 
  static void dotProduct( localIndex const size,
                          arrayView1d< real32 > const & vector1,
                          arrayView1d< real32 > const & vector2,
                          real64 & res )
  {
 
    RAJA::ReduceSum< parallelDeviceReduce, real64 > tmp( 0.0 );
    forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a )
    { tmp += vector1[a] * vector2[a]; } );
 
    res = tmp.get();
  }
 
  GEOS_HOST_DEVICE
  static R1Tensor computeDASVector( real64 const dip, real64 const azimuth )
  {
    real64 cd = cos( dip );
    real64 v1 = cd * cos( azimuth );
    real64 v2 = cd * sin( azimuth );
    real64 v3 = sin( dip );
    R1Tensor dasVector = {v1, v2, v3};
    return dasVector;
  }
 
  template< typename REAL >
  GEOS_HOST_DEVICE static REAL
  computeReferenceLengthForPenalty( arraySlice1d< localIndex const, cells::NODE_MAP_USD - 1 > const elemsToNodes,
                                    arrayView2d< REAL const, nodes::REFERENCE_POSITION_USD > const nodeCoords )
  {
    // Loop over the element faces
    REAL hs = 0;
    REAL spxf[3][2]{};
    REAL v1[3]{};
    REAL v2[3]{};
    REAL spx[3][3]{};
    for( int i = 0; i < 4; i++ )
    {
      REAL cross[3]{};
      int cf = 0;
      int jstart = (i + 1) % 4;
      int j = (jstart + 1) % 4;
      do
      {
        for( int k = 0; k < 3; k++ )
        {
          spxf[k][cf] = nodeCoords( elemsToNodes[j], k ) - nodeCoords( elemsToNodes[jstart], k );
        }
        j = (j + 1) % 4;
        cf++;
      }
      while (i != j);
 
      v1[0] = spxf[0][0];
      v1[1] = spxf[1][0];
      v1[2] = spxf[2][0];
      v2[0] = spxf[0][1];
      v2[1] = spxf[1][1];
      v2[2] = spxf[2][1];
      LvArray::tensorOps::crossProduct( cross, v1, v2 );
      hs = hs + LvArray::math::sqrt( LvArray::tensorOps::l2NormSquared< 3 >( cross ));
    }
 
    for( int i = 0; i < 3; i++ )
    {
      for( int k = 0; k < 3; k++ )
      {
        spx[i][k] = nodeCoords( elemsToNodes[i], k ) - nodeCoords( elemsToNodes[3], k );
      }
    }
 
 
    return LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( spx ) / hs );
  }
};
 
ENUM_STRINGS( WaveSolverUtils::DASType,
              "none",
              "dipole",
              "strainIntegration" );
 
ENUM_STRINGS( WaveSolverUtils::AttenuationType,
              "none",
              "sls" );
 
} /* namespace geos */
 
#endif /* GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERUTILS_HPP_ */