feat: MGR strategy for single phase flow + ALM poromechanics

src/coreComponents/common/format/Format.hpp

@@ -18,6 +18,8 @@
 
 #include <type_traits>
 #include <optional>
+#include <sstream>
+#include <string>
 
 #define GEOS_USE_FMT
 

src/coreComponents/linearAlgebra/CMakeLists.txt

@@ -73,6 +73,7 @@ set( linearAlgebra_headers
      utilities/ComponentMask.hpp
      utilities/InverseNormalOperator.hpp
      utilities/LAIHelperFunctions.hpp
+     utilities/SparsityPatternUtilities.hpp
      utilities/LinearSolverParameters.hpp
      utilities/LinearSolverResult.hpp
      utilities/NormalOperator.hpp
@@ -168,6 +169,8 @@ if( ENABLE_HYPRE )
           interfaces/hypre/mgrStrategies/MultiphasePoromechanics.hpp
           interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp
           interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp
+          interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp
+          interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp
           interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp
           interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp
           interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp

src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp

@@ -93,6 +93,15 @@ class LinearSolverBase : public PreconditionerBase< LAI >
     GEOS_UNUSED_VAR( context );
   }
 
+  /**
+   * @brief Supply near-null-space vectors used to construct physics-aware preconditioners.
+   * @param nearNullKernel Distributed vectors spanning the near null space.
+   */
+  virtual void setNearNullKernel( arrayView1d< Vector const > const & nearNullKernel )
+  {
+    GEOS_UNUSED_VAR( nearNullKernel );
+  }
+
   /**
    * @brief @return parameters of the solver.
    */

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

@@ -37,6 +37,8 @@
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanics.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp"
+#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp"
+#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp"
@@ -188,6 +190,16 @@ void hypre::mgr::createMGR( LinearSolverParameters const & params,
       setStrategy< SinglePhasePoromechanicsConformingFractures >( params.mgr, numComponentsPerField, precond, mgrData );
       break;
     }
+    case LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM:
+    {
+      setStrategy< SinglePhasePoromechanicsConformingFracturesALM >( params.mgr, numComponentsPerField, precond, mgrData );
+      break;
+    }
+    case LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALMReservoirFVM:
+    {
+      setStrategy< SinglePhasePoromechanicsConformingFracturesALMReservoirFVM >( params.mgr, numComponentsPerField, precond, mgrData );
+      break;
+    }
     case LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsReservoirFVM:
     {
       setStrategy< SinglePhasePoromechanicsReservoirFVM >( params.mgr, numComponentsPerField, precond, mgrData );

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

@@ -39,6 +39,7 @@ struct HypreMGRData
   array1d< HYPRE_Int > pointMarkers;  ///< array1d of unique tags for local degrees of freedom
   HyprePrecWrapper coarseSolver;      ///< MGR coarse solver pointer and functions
   HyprePrecWrapper mechSolver;        ///< MGR mechanics fine solver pointer and functions
+  HyprePrecWrapper nestedSolver;      ///< Optional nested MGR F-relaxation wrapper
 };
 
 namespace hypre
@@ -98,6 +99,21 @@ class MGRStrategyBase
   // HYPRE_Int m_numRestrictSweeps{ -1 }; ///< Number of restrict sweeps
   // HYPRE_Int m_numInterpSweeps{ -1 };   ///< Number of interpolation sweeps
 
+  /**
+   * @brief Total number of dof labels, i.e. the sum of all fields' components.
+   * @param numComponentsPerField number of components of each field
+   * @return the total number of blocks
+   */
+  static HYPRE_Int totalNumBlocks( arrayView1d< int const > const & numComponentsPerField )
+  {
+    HYPRE_Int result = 0;
+    for( localIndex i = 0; i < numComponentsPerField.size(); ++i )
+    {
+      result += LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[i] );
+    }
+    return result;
+  }
+
   /**
    * @brief Constructor.
    * @param numBlocks number of blocks
@@ -196,6 +212,58 @@ class MGRStrategyBase
     solver.destroy = HYPRE_BoomerAMGDestroy;
   }
 
+  /**
+   * @brief Set up one of the two BoomerAMG instances used by the fully
+   *        coupled single-phase ALM hierarchy.
+   * @param solver solver wrapper to initialize
+   * @param separateComponents whether displacement components are filtered
+   * @param bubbleCoarse true for the inner (bubble-displacement) coarse solve
+   *
+   * The CPU values mirror the reference nested-MGR YAML.  The device branch
+   * retains the device-safe smoother choices used by the existing strategies.
+   */
+  void setALMDisplacementAMG( HyprePrecWrapper & solver,
+                              integer const separateComponents,
+                              bool const bubbleCoarse )
+  {
+    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_BoomerAMGSetPrintLevel( solver.ptr, 0 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( solver.ptr, 1.0 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( solver.ptr, bubbleCoarse ? 0.75 : 0.8 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( solver.ptr, 3 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetFilterFunctions( solver.ptr, bubbleCoarse ? 0 : separateComponents ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPMaxElmts( solver.ptr, bubbleCoarse ? 10 : 20 ) );
+
+    if( !bubbleCoarse )
+    {
+      GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggNumLevels( solver.ptr, 1 ) );
+      GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( solver.ptr,
+                                                           hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::Falgout ) ) );
+    }
+
+#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 ) ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( solver.ptr, bubbleCoarse ? 1 : 2 ) );
+#else
+    HYPRE_Int constexpr l1SymmetricHybridGaussSeidel = 89;
+    HYPRE_Int constexpr gaussianElimination = 9;
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, l1SymmetricHybridGaussSeidel, 1 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, l1SymmetricHybridGaussSeidel, 2 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, gaussianElimination, 3 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( solver.ptr, bubbleCoarse ? 1 : 2 ) );
+    GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( solver.ptr, 0 ) );
+#endif
+
+    solver.setup = HYPRE_BoomerAMGSetup;
+    solver.solve = HYPRE_BoomerAMGSolve;
+    solver.destroy = HYPRE_BoomerAMGDestroy;
+  }
+
   /**
    * @brief Set up BoomerAMG to perform the solve for the pressure system
    * @param solver the solver wrapper

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

@@ -453,6 +453,10 @@ void HyprePreconditioner::clear()
   {
     GEOS_LAI_CHECK_ERROR( m_mgrData->mechSolver.destroy( m_mgrData->mechSolver.ptr ) );
   }
+  if( m_mgrData && m_mgrData->nestedSolver.ptr && m_mgrData->nestedSolver.destroy )
+  {
+    GEOS_LAI_CHECK_ERROR( m_mgrData->nestedSolver.destroy( m_mgrData->nestedSolver.ptr ) );
+  }
   m_precond.reset();
   m_mgrData.reset();
 }