PhysicsSolverBaseKernels.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 Total, S.A
 * 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_SOLVERBASEKERNELS_HPP
#define GEOS_PHYSICSSOLVERS_SOLVERBASEKERNELS_HPP
 
#include "common/format/EnumStrings.hpp"
#include "common/DataTypes.hpp"
#include "common/MpiWrapper.hpp"
 
namespace geos
{
 
namespace physicsSolverBaseKernels
{
 
/******************************** ResidualNormKernelBase ********************************/
 
template< integer NUM_NORM >
class ResidualNormKernelBase
{
public:
 
 
  static constexpr integer numNorm = NUM_NORM;
 
  ResidualNormKernelBase( globalIndex const rankOffset,
                          arrayView1d< real64 const > const & localResidual,
                          arrayView1d< globalIndex const > const & dofNumber,
                          arrayView1d< localIndex const > const & ghostRank,
                          real64 const minNormalizer ):
    m_rankOffset( rankOffset ),
    m_localResidual( localResidual ),
    m_dofNumber( dofNumber ),
    m_ghostRank( ghostRank ),
    m_minNormalizer( minNormalizer )
  {}
 
  struct LinfStackVariables
  {
    localIndex localRow;
 
    real64 localValue[numNorm]{};
  };
 
  struct L2StackVariables : public LinfStackVariables
  {
    real64 localNormalizer[numNorm]{};
  };
 
 
  GEOS_HOST_DEVICE
  integer ghostRank( localIndex const i ) const
  { return m_ghostRank( i ); }
 
  GEOS_HOST_DEVICE
  virtual void setupLinf( localIndex const i,
                          LinfStackVariables & stack ) const
  {
    stack.localRow = m_dofNumber[i] - m_rankOffset;
  }
 
  GEOS_HOST_DEVICE
  virtual void setupL2( localIndex const i,
                        L2StackVariables & stack ) const
  {
    stack.localRow = m_dofNumber[i] - m_rankOffset;
  }
 
 
  GEOS_HOST_DEVICE
  virtual void computeLinf( localIndex const i,
                            LinfStackVariables & stack ) const = 0;
 
  GEOS_HOST_DEVICE
  virtual void computeL2( localIndex const i,
                          L2StackVariables & stack ) const = 0;
 
  template< typename POLICY, typename KERNEL_TYPE >
  static void
  launchLinf( localIndex const size,
              KERNEL_TYPE const & kernelComponent,
              real64 (& residualNorm)[numNorm] )
  {
    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > localResidualNorm[numNorm]{};
 
    forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const i )
    {
      if( kernelComponent.ghostRank( i ) >= 0 )
      {
        return;
      }
 
      typename KERNEL_TYPE::LinfStackVariables stack;
      kernelComponent.setupLinf( i, stack );
      kernelComponent.computeLinf( i, stack );
 
      for( integer j = 0; j < numNorm; ++j )
      {
        localResidualNorm[j].max( LvArray::math::abs( stack.localValue[j] ) );
      }
    } );
 
    for( integer j = 0; j < numNorm; ++j )
    {
      residualNorm[j] = localResidualNorm[j].get();
    }
  }
 
  template< typename POLICY, typename KERNEL_TYPE >
  static void
  launchL2( localIndex const size,
            KERNEL_TYPE const & kernelComponent,
            real64 (& residualNorm)[numNorm],
            real64 (& residualNormalizer)[numNorm] )
  {
    RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > localResidualNorm[numNorm]{};
    RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > localResidualNormalizer[numNorm]{};
 
    forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const i )
    {
      if( kernelComponent.ghostRank( i ) >= 0 )
      {
        return;
      }
 
      typename KERNEL_TYPE::L2StackVariables stack;
      kernelComponent.setupL2( i, stack );
      kernelComponent.computeL2( i, stack );
 
      for( integer j = 0; j < numNorm; ++j )
      {
        localResidualNorm[j] += stack.localValue[j];
        localResidualNormalizer[j] += stack.localNormalizer[j];
      }
    } );
 
    for( integer j = 0; j < numNorm; ++j )
    {
      residualNorm[j] = localResidualNorm[j].get();
      residualNormalizer[j] = localResidualNormalizer[j].get();
    }
  }
 
 
protected:
 
  globalIndex const m_rankOffset;
 
  arrayView1d< real64 const > const m_localResidual;
 
  arrayView1d< globalIndex const > const m_dofNumber;
 
  arrayView1d< integer const > const m_ghostRank;
 
  real64 const m_minNormalizer;
 
};
 
class LinfResidualNormHelper
{
public:
 
  template< integer NUM_NORM >
  static void updateLocalNorm( real64 const (&subRegionResidualNorm)[NUM_NORM],
                               array1d< real64 > & localResidualNorm )
  {
    for( integer i = 0; i < NUM_NORM; ++i )
    {
      if( subRegionResidualNorm[i] > localResidualNorm[i] )
      {
        localResidualNorm[i] = subRegionResidualNorm[i];
      }
    }
  }
 
  static void computeGlobalNorm( real64 const & localResidualNorm,
                                 real64 & globalResidualNorm )
  {
    globalResidualNorm = MpiWrapper::max( localResidualNorm );
  }
 
  static void computeGlobalNorm( array1d< real64 > const & localResidualNorm,
                                 array1d< real64 > & globalResidualNorm )
  {
    MpiWrapper::allReduce( localResidualNorm.data(),
                           globalResidualNorm.data(),
                           localResidualNorm.size(),
                           MpiWrapper::getMpiOp( MpiWrapper::Reduction::Max ),
                           MPI_COMM_GEOS );
  }
};
 
class L2ResidualNormHelper
{
public:
 
  template< integer NUM_NORM >
  static void updateLocalNorm( real64 const (&subRegionResidualNorm)[NUM_NORM],
                               real64 const (&subRegionResidualNormalizer)[NUM_NORM],
                               array1d< real64 > & localResidualNorm,
                               array1d< real64 > & localResidualNormalizer )
  {
    for( integer i = 0; i < NUM_NORM; ++i )
    {
      localResidualNorm[i] += subRegionResidualNorm[i];
      localResidualNormalizer[i] += subRegionResidualNormalizer[i];
    }
  }
 
  static void computeGlobalNorm( real64 const & localResidualNorm,
                                 real64 const & localResidualNormalizer,
                                 real64 & globalResidualNorm )
  {
    globalResidualNorm = sqrt( MpiWrapper::sum( localResidualNorm ) ) / sqrt( MpiWrapper::sum( localResidualNormalizer ) );
  }
 
  static void computeGlobalNorm( array1d< real64 > const & localResidualNorm,
                                 array1d< real64 > const & localResidualNormalizer,
                                 array1d< real64 > & globalResidualNorm )
  {
    array1d< real64 > sumLocalResidualNorm( localResidualNorm.size() );
    array1d< real64 > sumLocalResidualNormalizer( localResidualNormalizer.size() );
    MpiWrapper::allReduce( localResidualNorm.data(),
                           sumLocalResidualNorm.data(),
                           localResidualNorm.size(),
                           MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
                           MPI_COMM_GEOS );
    MpiWrapper::allReduce( localResidualNormalizer.data(),
                           sumLocalResidualNormalizer.data(),
                           localResidualNormalizer.size(),
                           MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
                           MPI_COMM_GEOS );
    for( integer i = 0; i < localResidualNorm.size(); ++i )
    {
      globalResidualNorm[i] = sqrt( sumLocalResidualNorm[i] ) / sqrt( sumLocalResidualNormalizer[i] );
    }
  }
 
};
 
enum class NormType : integer
{
  Linf,  
  L2     
};
 
ENUM_STRINGS( NormType,
              "Linfinity",
              "L2" );
 
 
} // namespace physicsSolverBaseKernels
 
} // namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_SOLVERBASEKERNELS_HPP