WaveSolverBase.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_WAVESOLVERBASE_HPP_
#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERBASE_HPP_
 
 
#include "mesh/MeshFields.hpp"
#include "physicsSolvers/PhysicsSolverBase.hpp"
#include "common/LifoStorage.hpp"
#include "functions/TableFunction.hpp"
#if !defined( GEOS_USE_HIP )
#include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp"
#endif
#include "mesh/DomainPartition.hpp"
#include "WaveSolverUtils.hpp"
 
#if !defined( GEOS_USE_HIP )
#define SEM_FE_TYPES \
  finiteElement::Q1_Hexahedron_Lagrange_GaussLobatto, \
  finiteElement::Q2_Hexahedron_Lagrange_GaussLobatto, \
  finiteElement::Q3_Hexahedron_Lagrange_GaussLobatto, \
  finiteElement::Q4_Hexahedron_Lagrange_GaussLobatto, \
  finiteElement::Q5_Hexahedron_Lagrange_GaussLobatto
#else
#define SEM_FE_TYPES
#endif
 
#define SELECTED_FE_TYPES SEM_FE_TYPES
 
namespace geos
{
 
class WaveSolverBase : public PhysicsSolverBase
{
public:
 
  static constexpr real64 epsilonLoc = WaveSolverUtils::epsilonLoc;
  using EXEC_POLICY = WaveSolverUtils::EXEC_POLICY;
  using wsCoordType = WaveSolverUtils::wsCoordType;
 
  WaveSolverBase( const std::string & name,
                  Group * const parent );
 
  virtual ~WaveSolverBase() override;
 
  WaveSolverBase() = delete;
  WaveSolverBase( WaveSolverBase const & ) = delete;
  WaveSolverBase( WaveSolverBase && ) = default;
 
  WaveSolverBase & operator=( WaveSolverBase const & ) = delete;
  WaveSolverBase & operator=( WaveSolverBase && ) = delete;
 
  virtual void initializePreSubGroups() override;
 
  virtual real64 solverStep( real64 const & time_n,
                             real64 const & dt,
                             integer const cycleNumber,
                             DomainPartition & domain ) override;
 
 
  virtual real64 explicitStep( real64 const & time_n,
                               real64 const & dt,
                               integer const cycleNumber,
                               DomainPartition & domain ) override;
 
  struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct
  {
    static constexpr char const * sourceCoordinatesString() { return "sourceCoordinates"; }
    static constexpr char const * sourceValueString() { return "sourceValue"; }
 
    static constexpr char const * timeSourceFrequencyString() { return "timeSourceFrequency"; }
    static constexpr char const * timeSourceDelayString() { return "timeSourceDelay"; }
    static constexpr char const * rickerOrderString() { return "rickerOrder"; }
 
    static constexpr char const * receiverCoordinatesString() { return "receiverCoordinates"; }
 
    static constexpr char const * sourceNodeIdsString() { return "sourceNodeIds"; }
    static constexpr char const * sourceConstantsString() { return "sourceConstants"; }
    static constexpr char const * sourceIsAccessibleString() { return "sourceIsAccessible"; }
 
    static constexpr char const * receiverNodeIdsString() { return "receiverNodeIds"; }
    static constexpr char const * receiverConstantsString() {return "receiverConstants"; }
    static constexpr char const * receiverIsLocalString() { return "receiverIsLocal"; }
 
    static constexpr char const * outputSeismoTraceString() { return "outputSeismoTrace"; }
    static constexpr char const * dtSeismoTraceString() { return "dtSeismoTrace"; }
    static constexpr char const * indexSeismoTraceString() { return "indexSeismoTrace"; }
    static constexpr char const * forwardString() { return "forward"; }
    static constexpr char const * saveFieldsString() { return "saveFields"; }
    static constexpr char const * shotIndexString() { return "shotIndex"; }
    static constexpr char const * enableLifoString() { return "enableLifo"; }
    static constexpr char const * lifoSizeString() { return "lifoSize"; }
    static constexpr char const * lifoOnDeviceString() { return "lifoOnDevice"; }
    static constexpr char const * lifoOnHostString() { return "lifoOnHost"; }
 
    static constexpr char const * useDASString() { return "useDAS"; }
    static constexpr char const * linearDASSamplesString() { return "linearDASSamples"; }
    static constexpr char const * linearDASGeometryString() { return "linearDASGeometry"; }
    static constexpr char const * linearDASVectorXString() { return "linearDASVectorX"; }
    static constexpr char const * linearDASVectorYString() { return "linearDASVectorY"; }
    static constexpr char const * linearDASVectorZString() { return "linearDASVectorZ"; }
 
    static constexpr char const * usePMLString() { return "usePML"; }
    static constexpr char const * parametersPMLString() { return "parametersPML"; }
 
    static constexpr char const * receiverElemString() { return "receiverElem"; }
    static constexpr char const * receiverRegionString() { return "receiverRegion"; }
    static constexpr char const * freeSurfaceString() { return "FreeSurface"; }
 
    static constexpr char const * timestepStabilityLimitString() { return "timestepStabilityLimit"; }
    static constexpr char const * timeStepString() { return "timeStep"; }
 
    static constexpr char const * attenuationTypeString() { return "attenuationType"; }
    static constexpr char const * slsReferenceAngularFrequenciesString() { return "slsReferenceAngularFrequencies"; }
    static constexpr char const * slsAnelasticityCoefficientsString() { return "slsAnelasticityCoefficients"; }
 
    static constexpr char const * sourceWaveletTableNames() { return "sourceWaveletTableNames"; }
  };
 
  void reinit() override final;
 
  SortedArray< localIndex > const & getSolverNodesSet() { return m_solverTargetNodesSet; }
 
  void computeTargetNodeSet( arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes,
                             localIndex const subRegionSize,
                             localIndex const numQuadraturePointsPerElem );
 
  template< typename F, typename T >
  T computeGlobalMinOfElemField( DomainPartition & domain )
  {
    RAJA::ReduceMin< ReducePolicy< EXEC_POLICY >, T > minF( LvArray::NumericLimits< T >::max );
 
    forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                  MeshLevel & mesh,
                                                                  string_array const & regionNames )
    {
      mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
                                                                                            CellElementSubRegion & elementSubRegion )
      {
        arrayView1d< T const > const f = elementSubRegion.getField< F >();
        forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const e ) {
          minF.min( f[e] );
        } );
      } );
    } );
    T minFVal = minF.get();
    return MpiWrapper::min< T >( minFVal );
  }
 
  template< typename ... Qs >
  void initializeAnelasticityCoefficients( DomainPartition & domain )
  {
    // compute miniumum quality factor
    real32 minQVal = LvArray::NumericLimits< real32 >::max;
    ( (minQVal = std::min( minQVal, computeGlobalMinOfElemField< Qs, real32 >( domain ) ) ), ... );
    if( m_slsAnelasticityCoefficients.size( 0 ) == 1 && m_slsAnelasticityCoefficients[ 0 ] < 0 )
    {
      m_slsAnelasticityCoefficients[ 0 ] = 2.0 * minQVal / ( minQVal - 1.0 );
    }
    // test if anelasticity is too high and artifacts could appear
    real32 ySum = 0.0;
    for( integer l = 0; l < m_slsAnelasticityCoefficients.size( 0 ); l++ )
    {
      ySum += m_slsAnelasticityCoefficients[ l ];
    }
    GEOS_WARNING_IF( ySum > minQVal, "The anelasticity parameters are too high for the given quality factor. This could lead to solution artifacts such as zero-velocity waves." );
  }
 
protected:
 
  virtual void postInputInitialization() override;
 
  bool directoryExists( std::string const & directoryName );
 
  virtual void applyFreeSurfaceBC( real64 const time, DomainPartition & domain ) = 0;
 
  virtual real64 computeTimeStep( real64 & dtOut ) = 0;
 
  virtual void initializePML() = 0;
 
  virtual void incrementIndexSeismoTrace( real64 const time_n );
 
  virtual void computeAllSeismoTraces( real64 const time_n,
                                       real64 const dt,
                                       arrayView1d< real32 const > const var_np1,
                                       arrayView1d< real32 const > const var_n,
                                       arrayView2d< real32 > varAtReceivers,
                                       arrayView1d< real32 > coeffs = {},
                                       bool add = false );
  virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex,
                                                 real64 const time_n,
                                                 real64 const dt,
                                                 arrayView2d< real32 const > const var_np1,
                                                 arrayView2d< real32 const > const var_n,
                                                 arrayView2d< real32 > varAtReceivers );
 
  virtual void applyPML( real64 const time, DomainPartition & domain ) = 0;
 
  virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, string_array const & regionNames ) = 0;
 
  virtual real64 explicitStepForward( real64 const & time_n,
                                      real64 const & dt,
                                      integer const cycleNumber,
                                      DomainPartition & domain,
                                      integer const computeGradient ) = 0;
  virtual real64 explicitStepBackward( real64 const & time_n,
                                       real64 const & dt,
                                       integer const cycleNumber,
                                       DomainPartition & domain,
                                       integer const computeGradient ) = 0;
 
 
  virtual void registerDataOnMesh( Group & meshBodies ) override;
 
  localIndex getNumNodesPerElem();
 
  array2d< real32 > m_sourceValue;
 
  array2d< real64 > m_sourceCoordinates;
 
  real32 m_timeSourceFrequency;
 
  real32 m_timeSourceDelay;
 
  array2d< real64 > m_receiverCoordinates;
 
  localIndex m_rickerOrder;
 
  integer m_outputSeismoTrace;
 
  real64 m_dtSeismoTrace;
 
  localIndex m_indexSeismoTrace;
 
  localIndex m_nsamplesSeismoTrace;
 
  WaveSolverUtils::DASType m_useDAS;
 
  integer m_linearDASSamples;
 
  array2d< real64 > m_linearDASGeometry;
 
  array1d< real32 > m_linearDASVectorX;
 
  array1d< real32 > m_linearDASVectorY;
 
  array1d< real32 > m_linearDASVectorZ;
 
  WaveSolverUtils::AttenuationType m_attenuationType;
 
  array1d< real32 > m_slsReferenceAngularFrequencies;
 
  array1d< real32 > m_slsAnelasticityCoefficients;
 
  localIndex m_forward;
 
  integer m_saveFields;
 
  // Indicate the current shot computed for naming saved temporary data
  integer m_shotIndex;
 
  integer m_usePML;
 
  integer m_timestepStabilityLimit;
 
  //Time step computed with power iteration
  real64 m_timeStep;
 
  array2d< localIndex > m_sourceNodeIds;
 
  array2d< real64 > m_sourceConstants;
 
  array1d< localIndex > m_sourceIsAccessible;
 
  array2d< localIndex > m_receiverNodeIds;
 
  array2d< real64 > m_receiverConstants;
 
  array1d< localIndex > m_receiverIsLocal;
 
  array1d< localIndex > m_receiverElem;
 
  array1d< localIndex > m_receiverRegion;
 
  localIndex m_enableLifo;
 
  localIndex m_lifoSize;
 
  localIndex m_lifoOnDevice;
 
  localIndex m_lifoOnHost;
 
  std::unique_ptr< LifoStorage< real32, localIndex > > m_lifo;
 
  SortedArray< localIndex > m_solverTargetNodesSet;
 
  string_array m_sourceWaveletTableNames;
 
  bool m_useSourceWaveletTables;
 
  array1d< TableFunction::KernelWrapper > m_sourceWaveletTableWrappers;
 
  struct parametersPML
  {
    R1Tensor32 xMinPML;
 
    R1Tensor32 xMaxPML;
 
    real32 reflectivityPML;
 
    R1Tensor32 thicknessMinXYZPML;
    R1Tensor32 thicknessMaxXYZPML;
 
    R1Tensor32 waveSpeedMinXYZPML;
    R1Tensor32 waveSpeedMaxXYZPML;
  };
 
};
 
namespace fields
{
using reference32Type = array2d< WaveSolverUtils::wsCoordType, nodes::REFERENCE_POSITION_PERM >;
DECLARE_FIELD( referencePosition32,
               "referencePosition32",
               reference32Type,
               0,
               NOPLOT,
               WRITE_AND_READ,
               "Copy of the referencePosition from NodeManager in 32 bits integer" );
}
} /* namespace geos */
 
#endif /* GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERBASE_HPP_ */