feat: hypre improvements

src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.cpp

@@ -45,14 +45,23 @@ void HypreInterface::initialize()
                       "OMP_NUM_THREADS > 1 may not be optimal for certain hypre preconditioning options. " );
 #endif
 
-  HYPRE_Init();
+  HYPRE_Initialize();
 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
   HYPRE_SetExecutionPolicy( HYPRE_EXEC_DEVICE );
   HYPRE_SetSpGemmUseVendor( 0 );
   HYPRE_DeviceInitialize();
 #endif
   HYPRE_SetMemoryLocation( hypre::memoryLocation );
   HYPRE_SetPrintErrorMode( 1 );
+
+#if defined(HYPRE_USING_UMPIRE)
+  HYPRE_SetUmpireUMPoolName( "HYPRE_UM" );
+  HYPRE_SetUmpireHostPoolName( "HYPRE_HOST" );
+  HYPRE_SetUmpireDevicePoolName( "HYPRE_DEVICE" );
+  HYPRE_SetUmpirePinnedPoolName( "HYPRE_PINNED" );
+#endif
+
+  HYPRE_SetLogLevel( getenv( "HYPRE_LOG_LEVEL" ) ? atoi( getenv( "HYPRE_LOG_LEVEL" ) ) : 0 );
 }
 
 void HypreInterface::finalize()

src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.hpp

@@ -91,11 +91,7 @@ class MGRStrategyBase
   MGRCoarseGridMethod m_levelCoarseGridMethod[numLevels];     ///< Coarse grid method for each level
   MGRGlobalSmootherType m_levelGlobalSmootherType[numLevels]; ///< Global smoother type for each level
   HYPRE_Int m_levelGlobalSmootherIters[numLevels]{ -1 };      ///< Number of global smoother iterations for each level
-#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CPU
   HYPRE_Real m_coarseGridThreshold{ 1.0e-20 };                ///< Coarse grid truncation threshold
-#else
-  HYPRE_Real m_coarseGridThreshold{ 0.0 };                    ///< Coarse grid truncation threshold
-#endif
 
   // TODO: the following options are currently commented out in MGR's code.
   //       Let's consider their use when re-enable in hypre
@@ -167,21 +163,25 @@ class MGRStrategyBase
     GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetLevelSmoothIters( precond.ptr, m_levelGlobalSmootherIters ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetTruncateCoarseGridThreshold( precond.ptr, m_coarseGridThreshold ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetNonCpointsToFpoints( precond.ptr, 1 ));
+    GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetNonGalerkinMaxElmts( precond.ptr, 1 ));
   }
 
   /**
    * @brief Set up BoomerAMG to perform the solve for the displacement system
    * @param solver the solver wrapper
    */
-  void setDisplacementAMG( HyprePrecWrapper & solver )
+  void setDisplacementAMG( HyprePrecWrapper & solver,
+                           integer const & separateComponents )
   {
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &solver.ptr ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( solver.ptr, 0.0 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( solver.ptr, 1 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( solver.ptr, 1.0 ) );
-    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( solver.ptr, 0.8 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( solver.ptr, 0.6 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( solver.ptr, 0 ) );
-
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMinCoarseSize( solver.ptr, 512 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, hypre::getAMGCoarseType( LinearSolverParameters::AMG::CoarseType::gsElimWPivoting ), 3 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetInterpType(solver.ptr, hypre::getAMGInterpolationType( LinearSolverParameters::AMG::InterpType::modifiedExtendedE )) );
 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( solver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( solver.ptr, hypre::getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::chebyshev ) ) );
@@ -191,6 +191,7 @@ class MGRStrategyBase
 #endif
 
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( solver.ptr, 3 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetFilterFunctions( solver.ptr, separateComponents ) );
 
     solver.setup = HYPRE_BoomerAMGSetup;
     solver.solve = HYPRE_BoomerAMGSolve;
@@ -209,9 +210,11 @@ class MGRStrategyBase
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggNumLevels( solver.ptr, 1 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggPMaxElmts( solver.ptr, 16 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( solver.ptr, 0.0 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMinCoarseSize( solver.ptr, 1024 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, 99, 3 ) );
 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
-    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggNumLevels( solver.ptr, 0 ) );
-    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( solver.ptr, toUnderlying( AMGCoarseningType::PMIS ) ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggInterpType( solver.ptr, hypre::getAMGAggressiveInterpolationType( LinearSolverParameters::AMG::AggInterpType::modifiedExtendedE ) ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( solver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( solver.ptr, getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::l1jacobi ) ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( solver.ptr, 2 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( solver.ptr, 1.0 ) );
@@ -261,9 +264,10 @@ class MGRStrategyBase
    *       any other solver) for F-relaxation at any level
    */
   void setMechanicsFSolver( HyprePrecWrapper & precond,
-                            HypreMGRData & mgrData )
+                            HypreMGRData & mgrData,
+                            integer const & separateComponents )
   {
-    setDisplacementAMG( mgrData.mechSolver );
+    setDisplacementAMG( mgrData.mechSolver, separateComponents );
     HYPRE_MGRSetFSolver( precond.ptr, mgrData.mechSolver.solve, mgrData.mechSolver.setup, mgrData.mechSolver.ptr );
   }
 

src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp

@@ -340,32 +340,6 @@ void HyprePreconditioner::create( DofManager const * const dofManager )
 
 HypreMatrix const & HyprePreconditioner::setupPreconditioningMatrix( HypreMatrix const & mat )
 {
-  GEOS_MARK_FUNCTION;
-
-  if( m_params.preconditionerType == LinearSolverParameters::PreconditionerType::mgr && m_params.mgr.separateComponents )
-  {
-    GEOS_LAI_ASSERT_MSG( mat.dofManager() != nullptr, "MGR preconditioner requires a DofManager instance" );
-    HypreMatrix Pu;
-    HypreMatrix Auu;
-    {
-      Stopwatch timer( m_makeRestrictorTime );
-      mat.dofManager()->makeRestrictor( { { m_params.mgr.displacementFieldName, { 3, true } } }, mat.comm(), true, Pu );
-    }
-    {
-      Stopwatch timer( m_computeAuuTime );
-      mat.multiplyPtAP( Pu, Auu );
-    }
-    {
-      Stopwatch timer( m_componentFilterTime );
-      Auu.separateComponentFilter( m_precondMatrix, m_params.dofsPerNode );
-    }
-  }
-  else if( m_params.preconditionerType == LinearSolverParameters::PreconditionerType::amg && m_params.amg.separateComponents )
-  {
-    Stopwatch timer( m_componentFilterTime );
-    mat.separateComponentFilter( m_precondMatrix, m_params.dofsPerNode );
-    return m_precondMatrix;
-  }
   return mat;
 }
 
@@ -385,12 +359,6 @@ void HyprePreconditioner::setup( Matrix const & mat )
   {
     LvArray::system::FloatingPointExceptionGuard guard( FE_ALL_EXCEPT );
 
-    // Perform setup of the MGR mechanics F-solver with SDC matrix, if used
-    if( m_mgrData && m_mgrData->mechSolver.ptr && m_mgrData->mechSolver.setup )
-    {
-//      GEOS_LAI_CHECK_ERROR( m_mgrData->mechSolver.setup( m_mgrData->mechSolver.ptr, m_precondMatrix.unwrapped(), nullptr, nullptr ) );
-    }
-
     // Perform setup of the main solver, if needed
     if( m_precond->setup )
     {

src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.hpp

@@ -95,26 +95,6 @@ class HyprePreconditioner final : public PreconditionerBase< HypreInterface >
    */
   HyprePrecWrapper const & unwrapped() const;
 
-  /**
-   * @brief @return time spent setting up separate component matrix.
-   */
-  real64 componentFilterTime() const
-  {
-    return m_componentFilterTime;
-  }
-
-  /// @return time to construct restrictor matrix.
-  real64 makeRestrictorTime() const
-  {
-    return m_makeRestrictorTime;
-  }
-
-  /// @return time to apply restrictor matrix.
-  real64 computeAuuTime() const
-  {
-    return m_computeAuuTime;
-  }
-
 private:
 
   /**
@@ -144,15 +124,6 @@ class HyprePreconditioner final : public PreconditionerBase< HypreInterface >
 
   /// Null space vectors
   std::unique_ptr< HypreNullSpace > m_nullSpace;
-
-  /// Timing of separate component matrix construction
-  real64 m_componentFilterTime = 0.0;
-
-  /// Timing of the restrictor matrix construction
-  real64 m_makeRestrictorTime = 0.0;
-
-  /// Timing of the cost of applying the restrictor matrix to the system
-  real64 m_computeAuuTime = 0.0;
 };
 
 }

src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp

@@ -201,11 +201,7 @@ void HypreSolver::setup( HypreMatrix const & mat )
   clear();
   Base::setup( mat );
   Stopwatch timer( m_result.setupTime );
-
   m_precond.setup( mat );
-  m_componentFilterTime = m_precond.componentFilterTime();
-  m_makeRestrictorTime = m_precond.makeRestrictorTime();
-  m_computeAuuTime = m_precond.computeAuuTime();
 
   m_solver = std::make_unique< HypreSolverWrapper >();
   createHypreKrylovSolver( m_params, mat.comm(), *m_solver );
@@ -276,14 +272,17 @@ void HypreSolver::solve( HypreVector const & rhs,
 
   if( m_params.logLevel >= 1 )
   {
-    GEOS_LOG_RANK_0( "        Linear Solver | " << m_result.status <<
-                     " | Iterations: " << m_result.numIterations <<
-                     " | Final Rel Res: " << m_result.residualReduction <<
-                     " | Make Restrictor Time: " << m_makeRestrictorTime <<
-                     " | Compute Auu Time: " << m_computeAuuTime <<
-                     " | SC Filter Time: " << m_componentFilterTime <<
-                     " | Setup Time: " << m_result.setupTime << " s" <<
-                     " | Solve Time: " << m_result.solveTime << " s" );
+    HYPRE_BigInt global_num_rows, global_num_nonzeros;
+
+    // This involves an MPI collective call, and therefore we call it only when necessary
+    GEOS_LAI_CHECK_ERROR( HYPRE_IJMatrixGetGlobalInfo( matrix().unwrappedIJ(),
+                                                       &global_num_rows,
+                                                       &global_num_rows, // This is intentional and assuming the matrix is square
+                                                       &global_num_nonzeros) );
+
+    GEOS_LOG_RANK_0( GEOS_FMT( "        Linear Solver | {} | Unknowns: {} | Nonzeros: {} | Iterations: {} | Final Rel Res: {:.4e} | Setup Time: {:.3f} s | Solve Time: {:.3f} s",
+                     m_result.status, global_num_rows, global_num_nonzeros, m_result.numIterations,
+                     m_result.residualReduction, m_result.setupTime, m_result.solveTime ) );
   }
 }
 

src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.hpp

@@ -93,11 +93,6 @@ class HypreSolver final : public LinearSolverBase< HypreInterface >
 
   /// Pointers to hypre functions for the krylov solver
   std::unique_ptr< HypreSolverWrapper > m_solver;
-
-  /// Time of the most recent SC matrix construction
-  real64 m_componentFilterTime;
-  real64 m_makeRestrictorTime;
-  real64 m_computeAuuTime;
 };
 
 } // end geos namespace

src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp

@@ -399,6 +399,8 @@ inline HYPRE_Int getAMGCoarseType( LinearSolverParameters::AMG::CoarseType const
     { 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" );
 }

src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/HybridSinglePhasePoromechanics.hpp

@@ -84,10 +84,11 @@ class HybridSinglePhasePoromechanics : public MGRStrategyBase< 2 >
 
   /**
    * @brief Setup the MGR strategy.
+   * @param mgrParams MGR configuration parameters
    * @param precond preconditioner wrapper
    * @param mgrData auxiliary MGR data
    */
-  void setup( LinearSolverParameters::MGR const &,
+  void setup( LinearSolverParameters::MGR const & mgrParams,
               HyprePrecWrapper & precond,
               HypreMGRData & mgrData )
   {
@@ -96,7 +97,7 @@ class HybridSinglePhasePoromechanics : public MGRStrategyBase< 2 >
     GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetPMaxElmts( precond.ptr, 0 ));
 
     // Configure the BoomerAMG solver used as F-relaxation for the first level
-    setMechanicsFSolver( precond, mgrData );
+    setMechanicsFSolver( precond, mgrData, mgrParams.separateComponents );
 
     // Configure the BoomerAMG solver used as mgr coarse solver for the pressure reduced system
     setPressureAMG( mgrData.coarseSolver );

src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/Hydrofracture.hpp

@@ -70,17 +70,18 @@ class Hydrofracture : public MGRStrategyBase< 1 >
 
   /**
    * @brief Setup the MGR strategy.
+   * @param mgrParams MGR configuration parameters
    * @param precond preconditioner wrapper
    * @param mgrData auxiliary MGR data
    */
-  void setup( LinearSolverParameters::MGR const &,
+  void setup( LinearSolverParameters::MGR const & mgrParams,
               HyprePrecWrapper & precond,
               HypreMGRData & mgrData )
   {
     setReduction( precond, mgrData );
 
     // Configure the BoomerAMG solver used as F-relaxation for the first level
-    setMechanicsFSolver( precond, mgrData );
+    setMechanicsFSolver( precond, mgrData, mgrParams.separateComponents );
 
     // Configure the BoomerAMG solver used as mgr coarse solver for the pressure reduced system
     setPressureAMG( mgrData.coarseSolver );

src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/LagrangianContactMechanics.hpp

@@ -79,15 +79,15 @@ class LagrangianContactMechanics : public MGRStrategyBase< 1 >
    * @param precond preconditioner wrapper
    * @param mgrData auxiliary MGR data
    */
-  void setup( LinearSolverParameters::MGR const &,
+  void setup( LinearSolverParameters::MGR const & mgrParams,
               HyprePrecWrapper & precond,
               HypreMGRData & mgrData )
   {
     setReduction( precond, mgrData );
 
     // Configure the BoomerAMG solver used as mgr coarse solver for the displacement reduced system
     // (note that no separate displacement component approach is used here)
-    setDisplacementAMG( mgrData.coarseSolver );
+    setDisplacementAMG( mgrData.coarseSolver, mgrParams.separateComponents );
   }
 };
 

src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/MultiphasePoromechanics.hpp

@@ -72,20 +72,20 @@ class MultiphasePoromechanics : public MGRStrategyBase< 3 >
     setupLabels();
 
     // Level 0
-    m_levelFRelaxType[0]         = MGRFRelaxationType::amgVCycle;
-    m_levelFRelaxIters[0]        = 1;
-    m_levelInterpType[0]         = MGRInterpolationType::jacobi;
-    m_levelRestrictType[0]       = MGRRestrictionType::injection;
-    m_levelCoarseGridMethod[0]   = MGRCoarseGridMethod::nonGalerkin;
-    m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none;
+    m_levelFRelaxType[0]          = MGRFRelaxationType::amgVCycle;
+    m_levelFRelaxIters[0]         = 1;
+    m_levelInterpType[0]          = MGRInterpolationType::jacobi;
+    m_levelRestrictType[0]        = MGRRestrictionType::injection;
+    m_levelCoarseGridMethod[0]    = MGRCoarseGridMethod::nonGalerkin;
+    m_levelGlobalSmootherType[0]  = MGRGlobalSmootherType::none;
 
     // Level 1
-    m_levelFRelaxType[1]         = MGRFRelaxationType::jacobi;
-    m_levelFRelaxIters[1]        = 1;
-    m_levelInterpType[1]         = MGRInterpolationType::jacobi;
-    m_levelRestrictType[1]       = MGRRestrictionType::injection;
-    m_levelCoarseGridMethod[1]   = MGRCoarseGridMethod::galerkin;
-    m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::none;
+    m_levelFRelaxType[1]          = MGRFRelaxationType::jacobi;
+    m_levelFRelaxIters[1]         = 1;
+    m_levelInterpType[1]          = MGRInterpolationType::jacobi;
+    m_levelRestrictType[1]        = MGRRestrictionType::injection;
+    m_levelCoarseGridMethod[1]    = MGRCoarseGridMethod::galerkin;
+    m_levelGlobalSmootherType[1]  = MGRGlobalSmootherType::none;
 
     // Level 2
     m_levelFRelaxType[2]          = MGRFRelaxationType::none;
@@ -98,17 +98,18 @@ class MultiphasePoromechanics : public MGRStrategyBase< 3 >
 
   /**
    * @brief Setup the MGR strategy.
+   * @param mgrParams MGR configuration parameters
    * @param precond preconditioner wrapper
    * @param mgrData auxiliary MGR data
    */
-  void setup( LinearSolverParameters::MGR const &,
+  void setup( LinearSolverParameters::MGR const & mgrParams,
               HyprePrecWrapper & precond,
               HypreMGRData & mgrData )
   {
     setReduction( precond, mgrData );
 
     // Configure the BoomerAMG solver used as F-relaxation for the first level
-    setMechanicsFSolver( precond, mgrData );
+    setMechanicsFSolver( precond, mgrData, mgrParams.separateComponents );
 
     // Configure the BoomerAMG solver used as mgr coarse solver for the pressure reduced system
     setPressureAMG( mgrData.coarseSolver );

src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/MultiphasePoromechanicsReservoirFVM.hpp

@@ -118,7 +118,7 @@ class MultiphasePoromechanicsReservoirFVM : public MGRStrategyBase< 4 >
 
   /**
    * @brief Setup the MGR strategy.
-   * @param mgrParams parameters for the configuration of the MGR recipe
+   * @param mgrParams MGR configuration parameters
    * @param precond preconditioner wrapper
    * @param mgrData auxiliary MGR data
    */
@@ -136,7 +136,7 @@ class MultiphasePoromechanicsReservoirFVM : public MGRStrategyBase< 4 >
     setReduction( precond, mgrData );
 
     // Configure the BoomerAMG solver used as F-relaxation for the first level
-    setMechanicsFSolver( precond, mgrData );
+    setMechanicsFSolver( precond, mgrData, mgrParams.separateComponents );
 
     // Configure the BoomerAMG solver used as mgr coarse solver for the pressure reduced system
     setPressureAMG( mgrData.coarseSolver );