HypreUtils.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_LINEARALGEBRA_INTERFACES_HYPREUTILS_HPP_
#define GEOS_LINEARALGEBRA_INTERFACES_HYPREUTILS_HPP_
 
#include "common/DataTypes.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
#include "common/MpiWrapper.hpp"
 
#include "codingUtilities/Utilities.hpp"
#include "linearAlgebra/utilities/LinearSolverParameters.hpp"
 
#include <HYPRE_krylov.h>
#include <HYPRE_parcsr_ls.h>
 
#define GEOS_HYPRE_DEVICE GEOS_HOST_DEVICE
 
namespace geos
{
 
class HypreMatrix;
class HypreVector;
 
struct HyprePrecWrapper
{
  using SetupFunc = HYPRE_Int (*)( HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector );
 
  using SolveFunc = HYPRE_Int (*)( HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector );
 
  using DestroyFunc = HYPRE_Int (*)( HYPRE_Solver );
 
  HYPRE_Solver ptr{};    
  SetupFunc setup{};     
  SolveFunc solve{};     
  DestroyFunc destroy{}; 
};
 
namespace hypre
{
 
constexpr HYPRE_MemoryLocation getMemoryLocation( LvArray::MemorySpace const space )
{
  switch( space )
  {
    case hostMemorySpace: return HYPRE_MEMORY_HOST;
#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA
    case LvArray::MemorySpace::cuda: return HYPRE_MEMORY_DEVICE;
#endif
#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
    case  LvArray::MemorySpace::hip: return HYPRE_MEMORY_DEVICE;
#endif
    default: return HYPRE_MEMORY_HOST;
  }
}
 
constexpr LvArray::MemorySpace getLvArrayMemorySpace( HYPRE_MemoryLocation const location )
{
  switch( location )
  {
    case HYPRE_MEMORY_HOST: return hostMemorySpace;
#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA
    case HYPRE_MEMORY_DEVICE: return parallelDeviceMemorySpace;
#endif
#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
    case HYPRE_MEMORY_DEVICE: return parallelDeviceMemorySpace;
#endif
    default: return hostMemorySpace;
  }
}
 
#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
 
using execPolicy = parallelDevicePolicy<>;
constexpr LvArray::MemorySpace memorySpace = parallelDeviceMemorySpace;
constexpr HYPRE_MemoryLocation memoryLocation = HYPRE_MEMORY_DEVICE;
 
#else
 
using execPolicy = parallelHostPolicy;
constexpr LvArray::MemorySpace memorySpace = hostMemorySpace;
constexpr HYPRE_MemoryLocation memoryLocation = HYPRE_MEMORY_HOST;
 
#endif
 
// Check matching requirements on index/value types between GEOS and Hypre
 
// WARNING. We don't have consistent types between HYPRE_Int and localIndex.
//          Decision needs to be made either to use bigint option, or change
//          localIndex to int. We are getting away with this because we do not
//          pass ( localIndex * ) to hypre except when it is on the GPU, in
//          which case we are using int for localIndex.
#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
static_assert( sizeof( HYPRE_Int ) == sizeof( geos::localIndex ),
               "HYPRE_Int and geos::localIndex must have the same size" );
static_assert( std::is_signed< HYPRE_Int >::value == std::is_signed< geos::localIndex >::value,
               "HYPRE_Int and geos::localIndex must both be signed or unsigned" );
#endif
 
inline void checkDeviceErrors( char const * msg, char const * file, int const line )
{
#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA
  cudaError_t const err = cudaGetLastError();
  GEOS_ERROR_IF( err != cudaSuccess, GEOS_FMT( "Previous CUDA errors found: {} ({} at {}:{})", msg, cudaGetErrorString( err ), file, line ) );
#elif GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
  hipError_t const err = hipGetLastError();
  GEOS_UNUSED_VAR( msg, file, line ); // on crusher geos_error_if ultimately resolves to an assert, which drops the content on release
                                      // builds
  GEOS_ERROR_IF( err != hipSuccess, GEOS_FMT( "Previous HIP errors found: {} ({} at {}:{})", msg, hipGetErrorString( err ), file, line ) );
#else
  GEOS_UNUSED_VAR( msg, file, line );
#endif
}
 
#define GEOS_HYPRE_CHECK_DEVICE_ERRORS( msg ) ::geos::hypre::checkDeviceErrors( msg, __FILE__, __LINE__ )
 
static_assert( sizeof( HYPRE_BigInt ) == sizeof( geos::globalIndex ),
               "HYPRE_BigInt and geos::globalIndex must have the same size" );
 
static_assert( std::is_signed< HYPRE_BigInt >::value == std::is_signed< geos::globalIndex >::value,
               "HYPRE_BigInt and geos::globalIndex must both be signed or unsigned" );
 
static_assert( std::is_same< HYPRE_Real, geos::real64 >::value,
               "HYPRE_Real and geos::real64 must be the same type" );
 
inline HYPRE_BigInt * toHypreBigInt( geos::globalIndex * const index )
{
  return reinterpret_cast< HYPRE_BigInt * >(index);
}
 
inline HYPRE_BigInt const * toHypreBigInt( geos::globalIndex const * const index )
{
  return reinterpret_cast< HYPRE_BigInt const * >(index);
}
 
HYPRE_Vector parVectorToVectorAll( HYPRE_ParVector const vec );
 
HYPRE_Vector parVectorToVector( HYPRE_ParVector const vec, int const targetRank );
 
HYPRE_Int dummySetup( HYPRE_Solver,
                      HYPRE_ParCSRMatrix,
                      HYPRE_ParVector,
                      HYPRE_ParVector );
 
void fillKrylovDofLabels( HypreMatrix const & mat,
                          array1d< int > & labels );
 
void assignKrylovDofTags( arrayView1d< int const > const & labels,
                          MPI_Comm const comm,
                          array1d< HYPRE_Int > & tags,
                          HYPRE_Int & numTags );
 
namespace testing
{
 
void setKrylovDofLabels( arrayView1d< int const > const & labels );
 
void clearKrylovDofLabels();
 
}
 
HYPRE_Int SuperLUDistSolve( HYPRE_Solver solver,
                            HYPRE_ParCSRMatrix A,
                            HYPRE_ParVector b,
                            HYPRE_ParVector x );
 
HYPRE_Int SuperLUDistDestroy( HYPRE_Solver solver );
 
HYPRE_Int relaxationCreate( HYPRE_Solver & solver,
                            HYPRE_Int const type );
 
HYPRE_Int relaxationSetup( HYPRE_Solver solver,
                           HYPRE_ParCSRMatrix A,
                           HYPRE_ParVector b,
                           HYPRE_ParVector x );
 
HYPRE_Int relaxationSolve( HYPRE_Solver solver,
                           HYPRE_ParCSRMatrix A,
                           HYPRE_ParVector b,
                           HYPRE_ParVector x );
 
HYPRE_Int relaxationDestroy( HYPRE_Solver solver );
 
HYPRE_Int chebyshevCreate( HYPRE_Solver & solver,
                           HYPRE_Int const order,
                           HYPRE_Int const eigNumIter );
 
HYPRE_Int chebyshevSetup( HYPRE_Solver solver,
                          HYPRE_ParCSRMatrix A,
                          HYPRE_ParVector b,
                          HYPRE_ParVector x );
 
HYPRE_Int chebyshevSolve( HYPRE_Solver solver,
                          HYPRE_ParCSRMatrix A,
                          HYPRE_ParVector b,
                          HYPRE_ParVector x );
 
HYPRE_Int chebyshevDestroy( HYPRE_Solver solver );
 
inline HYPRE_Int getAMGCycleType( LinearSolverParameters::AMG::CycleType const & type )
{
  static map< LinearSolverParameters::AMG::CycleType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::AMG::CycleType::V, 1 },
    { LinearSolverParameters::AMG::CycleType::W, 2 },
  };
  return findOption( typeMap, type, "multigrid cycle", "HyprePreconditioner" );
}
 
inline HYPRE_Int getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType const & type )
{
  static map< LinearSolverParameters::AMG::SmootherType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::AMG::SmootherType::default_, -1 },
    { LinearSolverParameters::AMG::SmootherType::fgs, 3 },
    { LinearSolverParameters::AMG::SmootherType::bgs, 4 },
    { LinearSolverParameters::AMG::SmootherType::sgs, 6 },
    { LinearSolverParameters::AMG::SmootherType::jacobi, 7 },
    { LinearSolverParameters::AMG::SmootherType::l1sgs, 8 },
    { LinearSolverParameters::AMG::SmootherType::chebyshev, 16 },
    { LinearSolverParameters::AMG::SmootherType::l1jacobi, 18 },
  };
  return findOption( typeMap, type, "multigrid relaxation", "HyprePreconditioner" );
}
 
inline HYPRE_Int getAMGInterpolationType( LinearSolverParameters::AMG::InterpType const & type )
{
  static map< LinearSolverParameters::AMG::InterpType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::AMG::InterpType::default_, -1 },
    { LinearSolverParameters::AMG::InterpType::modifiedClassical, 0 },
    { LinearSolverParameters::AMG::InterpType::direct, 3 },
    { LinearSolverParameters::AMG::InterpType::multipass, 4 },
    { LinearSolverParameters::AMG::InterpType::extendedI, 6 },
    { LinearSolverParameters::AMG::InterpType::standard, 8 },
    { LinearSolverParameters::AMG::InterpType::extended, 14 },
    { LinearSolverParameters::AMG::InterpType::directBAMG, 15 },
    { LinearSolverParameters::AMG::InterpType::modifiedExtended, 16 },
    { LinearSolverParameters::AMG::InterpType::modifiedExtendedI, 17 },
    { LinearSolverParameters::AMG::InterpType::modifiedExtendedE, 18 },
  };
  return findOption( typeMap, type, "multigrid interpolation", "HyprePreconditioner" );
}
 
inline HYPRE_Int getAMGAggressiveInterpolationType( LinearSolverParameters::AMG::AggInterpType const & type )
{
  static map< LinearSolverParameters::AMG::AggInterpType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::AMG::AggInterpType::default_, -1 },
    { LinearSolverParameters::AMG::AggInterpType::extendedIStage2, 1 },
    { LinearSolverParameters::AMG::AggInterpType::standardStage2, 2 },
    { LinearSolverParameters::AMG::AggInterpType::extendedStage2, 3 },
    { LinearSolverParameters::AMG::AggInterpType::multipass, 4 },
    { LinearSolverParameters::AMG::AggInterpType::modifiedExtended, 5 },
    { LinearSolverParameters::AMG::AggInterpType::modifiedExtendedI, 6 },
    { LinearSolverParameters::AMG::AggInterpType::modifiedExtendedE, 7 },
    { LinearSolverParameters::AMG::AggInterpType::modifiedMultipass, 8 },
  };
  return findOption( typeMap, type, "multigrid aggressive interpolation", "HyprePreconditioner" );
}
 
inline HYPRE_Int getILUType( LinearSolverParameters::AMG::SmootherType const type )
{
  static map< LinearSolverParameters::AMG::SmootherType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::AMG::SmootherType::iluk, 0 },
    { LinearSolverParameters::AMG::SmootherType::ilut, 1 },
  };
  return findOption( typeMap, type, "ILU", "HyprePreconditioner" );
}
 
inline HYPRE_Int getAMGCoarseType( LinearSolverParameters::AMG::CoarseType const & type )
{
  static map< LinearSolverParameters::AMG::CoarseType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::AMG::CoarseType::default_, -1 },
    { LinearSolverParameters::AMG::CoarseType::fgs, 3 },
    { LinearSolverParameters::AMG::CoarseType::bgs, 4 },
    { LinearSolverParameters::AMG::CoarseType::sgs, 6 },
    { LinearSolverParameters::AMG::CoarseType::jacobi, 7 },
    { LinearSolverParameters::AMG::CoarseType::l1sgs, 8 },
    { LinearSolverParameters::AMG::CoarseType::direct, 9 },
    { LinearSolverParameters::AMG::CoarseType::chebyshev, 16 },
    { LinearSolverParameters::AMG::CoarseType::l1jacobi, 18 },
    { LinearSolverParameters::AMG::CoarseType::gsElimWPivoting, 99 },
    { LinearSolverParameters::AMG::CoarseType::gsElimWInverse, 199 },
  };
  return findOption( typeMap, type, "multigrid coarse solver", "HyprePreconditioner" );
}
 
inline HYPRE_Int getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType const & type )
{
  static map< LinearSolverParameters::AMG::CoarseningType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::AMG::CoarseningType::default_, -1 },
    { LinearSolverParameters::AMG::CoarseningType::CLJP, 0 },
    { LinearSolverParameters::AMG::CoarseningType::RugeStueben, 3 },
    { LinearSolverParameters::AMG::CoarseningType::Falgout, 6 },
    { LinearSolverParameters::AMG::CoarseningType::PMIS, 8 },
    { LinearSolverParameters::AMG::CoarseningType::HMIS, 10 },
  };
  return findOption( typeMap, type, "multigrid coarsening", "HyprePreconditioner" );
}
 
inline HYPRE_Int getRelaxationType( LinearSolverParameters::PreconditionerType const type )
{
  static map< LinearSolverParameters::PreconditionerType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::PreconditionerType::fgs, 3 },
    { LinearSolverParameters::PreconditionerType::bgs, 4 },
    { LinearSolverParameters::PreconditionerType::sgs, 6 },
    { LinearSolverParameters::PreconditionerType::jacobi, 7 },
    { LinearSolverParameters::PreconditionerType::l1sgs, 8 },
    { LinearSolverParameters::PreconditionerType::chebyshev, 16 },
    { LinearSolverParameters::PreconditionerType::l1jacobi, 18 },
  };
  return findOption( typeMap, type, "relaxation", "HyprePreconditioner" );
}
 
inline HYPRE_Int getILUType( LinearSolverParameters::PreconditionerType const type )
{
  static map< LinearSolverParameters::PreconditionerType, HYPRE_Int > const typeMap =
  {
    { LinearSolverParameters::PreconditionerType::iluk, 0 },
    { LinearSolverParameters::PreconditionerType::ilut, 1 },
  };
  return findOption( typeMap, type, "ILU", "HyprePreconditioner" );
}
 
enum class AMGCoarseningType : HYPRE_Int
{
  CLJP = 0,         
  Ruge_Stueben = 3, 
  Falgout = 6,      
  CLJPDebug = 7,    
  PMIS = 8,         
  PMISDebug = 9,    
  HMIS = 10,        
  CGC = 21,         
  CGC_E = 22        
};
 
enum class MGRInterpolationType : HYPRE_Int
{
  injection = 0,                      
  l1jacobi = 1,                       
  jacobi = 2,                         
  classicalModifiedInterpolation = 3, 
  approximateInverse = 4,             
  blockJacobi = 12                    
};
 
enum class MGRRestrictionType : HYPRE_Int
{
  injection = 0,           
  jacobi = 2,              
  approximateInverse = 3,  
  blockJacobi = 12,        
  cprLike = 13,            
  blockColLumped = 14,     
  partialColLumped = 15    
};
 
enum class MGRCoarseGridMethod : HYPRE_Int
{
  galerkin = 0,           
  nonGalerkin = 1,        
  cprLikeDiag = 2,        
  cprLikeBlockDiag = 3,   
  approximateInverse = 4  
};
 
enum class MGRFRelaxationType : HYPRE_Int
{
  none = -1,                        
  singleVCycleSmoother = 1,         
  amgVCycle = 2,                    
  forwardHybridGaussSeidel = 3,     
  backwardHybridGaussSeidel = 4,    
  hybridSymmetricGaussSeidel = 6,   
  jacobi = 7,                       
  l1hybridSymmetricGaussSeidel = 8, 
  gsElim = 9,                       
  l1forwardGaussSeidel = 13,        
  l1backwardGaussSeidel = 14,       
  l1jacobi = 18,                    
  ilu = 32,                         
  gsElimWPivoting = 99,             
  gsElimWInverse = 199              
};
 
enum class MGRGlobalSmootherType : HYPRE_Int
{
  none = -1,            
  blockJacobi = 0,      
  blockGaussSeidel = 1, 
  jacobi = 2,           
  ilu0 = 16             
};
 
} // namespace hypre
 
} // namespace geos
 
#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREUTILS_HPP_*/