feat: hypre improvements (#3339)

Co-authored-by: Matteo Cusini <49037133+CusiniM@users.noreply.github.com>

.devcontainer/devcontainer.json

@@ -2,7 +2,7 @@
     "build": {
         "dockerfile": "Dockerfile",
         "args": {
-            "GEOS_TPL_TAG": "284-535"
+            "GEOS_TPL_TAG": "286-547"
         }
     },
     "runArgs": [

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3434-8663-1be8ad2
+  baseline: integratedTests/baseline_integratedTests-pr3339-8707-7c55c70
 allow_fail:
   all: ''
   streak: ''

BASELINE_NOTES.md

@@ -6,6 +6,11 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+
+PR #3339 (2024-11-14)
+=====================
+Hypre improvements, rebaseline is due to field value change (amgNumFunctions).
+
 PR #3434 (2024-11-09)
 =====================
 Bugfix: Fixed output of ArrayOfArray objects to restart files.

src/coreComponents/common/format/StringUtilities.cpp

@@ -73,6 +73,29 @@ string removeStringAndFollowingContent( string_view const str,
   return string( newStr );
 }
 
+// Add comma separators for thousands
+template< typename T >
+string addCommaSeparators( T const & num )
+{
+  static_assert( std::is_integral< T >::value, "addCommaSeparators only supports integral types" );
+
+  string const numStr = std::to_string( num );
+  string result;
+
+  for( std::size_t i = 0; i < numStr.size(); ++i )
+  {
+    result += numStr[i];
+    if((numStr.size() - i - 1) % 3 == 0 && i != numStr.size() - 1 )
+    {
+      result += ",";
+    }
+  }
+  return result;
+}
+
+template string addCommaSeparators( localIndex const & num );
+template string addCommaSeparators( globalIndex const & num );
+
 // put definition here so we can control the allowable values of T and
 // modication of this function triggers a whole code recompile...which
 // should be avoided.

src/coreComponents/common/format/StringUtilities.hpp

@@ -227,6 +227,15 @@ string_view trimSpaces( string_view str );
 string removeStringAndFollowingContent( string_view str,
                                         string_view strToRemove );
 
+/**
+ * @brief Add comma separators to an integral number for readability.
+ * @tparam T the integral type of the number to format.
+ * @param[in] num the integral number to format.
+ * @return a string representation of the number with comma separators.
+ */
+template< typename T >
+string addCommaSeparators( T const & num );
+
 /**
  * @brief Take a string, and return a array1d with the cast values
  * @tparam T the type to which the string will be cast

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

@@ -40,19 +40,35 @@ namespace geos
 
 void HypreInterface::initialize()
 {
-#ifdef GEOS_USE_OPENMP
-  GEOS_LOG_RANK_0_IF( omp_get_max_threads()>1,
-                      "OMP_NUM_THREADS > 1 may not be optimal for certain hypre preconditioning options. " );
+#if defined(GEOS_USE_OPENMP) && defined(HYPRE_USING_OPENMP)
+  GEOS_LOG_RANK_0_IF( omp_get_max_threads() > 1,
+                      "\n"
+                      "********************************************************************\n"
+                      "*                                                                  *\n"
+                      "*    WARNING: OMP_NUM_THREADS > 1 MAY NOT BE OPTIMAL FOR CERTAIN   *\n"
+                      "*             HYPRE PRECONDITIONING OPTIONS!                       *\n"
+                      "*                                                                  *\n"
+                      "********************************************************************\n"
+                      );
 #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,24 @@ 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
+   * @param separateComponents flag controlling the use of the separate displacement component (SDC) approximation
    */
-  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_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 +190,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;
@@ -207,11 +207,12 @@ class MGRStrategyBase
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( solver.ptr, 0 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( solver.ptr, 1 ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggNumLevels( solver.ptr, 1 ) );
-    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggPMaxElmts( solver.ptr, 16 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggPMaxElmts( solver.ptr, 20 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggInterpType( solver.ptr, hypre::getAMGAggressiveInterpolationType( LinearSolverParameters::AMG::AggInterpType::multipass ) ) );
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( solver.ptr, 0.0 ) );
 #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 ) );
@@ -256,14 +257,16 @@ class MGRStrategyBase
    * @brief Set up BoomerAMG to perform the mechanics F-solve for the first F-relaxation
    * @param precond the preconditioner wrapper
    * @param mgrData auxiliary MGR data
+   * @param separateComponents flag controlling the use of the separate displacement component (SDC) approximation
    *
    * @note This function should be rethought once MGR allows for customizing boomerAMG (or
    *       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

@@ -88,7 +88,6 @@ void createAMG( LinearSolverParameters const & params,
   GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( precond.ptr, 0.0 ) );
   GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( precond.ptr, 1 ) );
   GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( precond.ptr, logLevel ) );
-  GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( precond.ptr, params.dofsPerNode ) );
 
   // Set maximum number of multigrid levels (default 25)
   GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxLevels( precond.ptr, LvArray::integerConversion< HYPRE_Int >( params.amg.maxLevels ) ) );
@@ -174,7 +173,9 @@ void createAMG( LinearSolverParameters const & params,
     GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPMaxElmts( precond.ptr, params.amg.interpolationMaxNonZeros ) );
   }
 
+  // Unknown-based AMG parameters
   GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( precond.ptr, params.amg.numFunctions ) );
+  GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetFilterFunctions( precond.ptr, params.amg.separateComponents ) );
 
   if( params.amg.aggressiveNumLevels )
   {
@@ -340,32 +341,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 +360,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,18 @@ 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, stringutilities::addCommaSeparators( global_num_rows ),
+                               stringutilities::addCommaSeparators( 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 );