feat: MGR reciepe for SinglePhase ALM

src/coreComponents/common/format/LogPart.hpp

@@ -197,8 +197,6 @@ void LogPart::addDescriptionBySection( Description & description, FormattedDescr
   size_t & formattedDescriptionMaxWidth = formattedDescription.m_maxValueWidth;
   size_t & formattedDescriptionNameWidth = formattedDescription.m_maxNameWidth;
   ( [&] {
-    static_assert( has_formatter_v< decltype(args) >,
-                   "Argument passed cannot be converted to string" );
     string const value = GEOS_FMT( "{}", args );
 
     stdVector< string_view > dividedDescriptionValues =

src/coreComponents/common/format/table/TableData.hpp

@@ -267,7 +267,6 @@ void TableData::addRow( Args const &... args )
 {
   stdVector< CellData > cells;
   ( [&] {
-    static_assert( has_formatter_v< decltype(args) > || isCellType< std::decay_t< decltype(args) > >, "Argument passed in addRow cannot be converted to string nor a CellType" );
     if constexpr (std::is_same_v< Args, CellType >) {
       cells.push_back( { args, string() } );
     }
@@ -297,7 +296,6 @@ void TableData::addRow( Args const &... args )
 template< typename T >
 void TableData2D::addCell( real64 const rowValue, real64 const columnValue, T const & value )
 {
-  static_assert( has_formatter_v< decltype(value) >, "Argument passed in addCell cannot be converted to string" );
   m_columnValues.insert( columnValue );
   m_data.get_inserted( rowValue ).get_inserted( columnValue ) =  GEOS_FMT( "{}", value );
 }

src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp

@@ -32,6 +32,7 @@
 #include <vtkSmartPointer.h>
 #include <vtkThreshold.h>
 #include <vtkUnstructuredGrid.h>
+#include <vtkVersionMacros.h>
 #include <vtkXMLUnstructuredGridWriter.h>
 #include <vtkAggregateDataSetFilter.h>
 // System includes
@@ -308,7 +309,9 @@ getWell( WellElementSubRegion const & subRegion,
   localIndex const numPoints = subRegion.size() > 0 ? subRegion.size() + 1 : 0;
   points->SetNumberOfPoints( numPoints );
   auto cellsArray = vtkSmartPointer< vtkCellArray >::New();
+#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 )
   cellsArray->SetNumberOfCells( subRegion.size() );
+#endif
   localIndex const numberOfNodesPerElement = subRegion.numNodesPerElement();
   GEOS_ERROR_IF_NE( numberOfNodesPerElement, 2 );
   stdVector< vtkIdType > connectivity( numberOfNodesPerElement );
@@ -358,7 +361,9 @@ getSurface( FaceElementSubRegion const & subRegion,
   auto & faceToNodes = faceManager.nodeList();
 
   auto cellArray = vtkSmartPointer< vtkCellArray >::New();
+#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 )
   cellArray->SetNumberOfCells( subRegion.size() );
+#endif
   stdVector< int > cellTypes;
   cellTypes.reserve( subRegion.size() );
 
@@ -599,7 +604,9 @@ static ParticleData
 getVtkCells( ParticleRegion const & region )
 {
   vtkSmartPointer< vtkCellArray > cellsArray = vtkCellArray::New();
+#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 )
   cellsArray->SetNumberOfCells( region.getNumberOfParticles< ParticleRegion >() );
+#endif
   stdVector< int > cellType;
   cellType.reserve( region.getNumberOfParticles< ParticleRegion >() );
 

src/coreComponents/linearAlgebra/CMakeLists.txt

@@ -168,6 +168,7 @@ 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/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,7 @@
 #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/SinglePhasePoromechanicsReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp"
@@ -188,6 +189,11 @@ 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::singlePhasePoromechanicsReservoirFVM:
     {
       setStrategy< SinglePhasePoromechanicsReservoirFVM >( params.mgr, numComponentsPerField, precond, mgrData );

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

@@ -23,6 +23,7 @@
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseHybridFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanics.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp"
+#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp"
@@ -40,6 +41,7 @@
 
 #include <algorithm>
 #include <cctype>
+#include <cstdio>
 #include <cstdlib>
 #include <fstream>
 #include <numeric>
@@ -57,12 +59,26 @@ namespace
 {
 
 void checkHypredriveCall( uint32_t const errorCode,
-                          char const * const call )
+                          char const * const call,
+                          std::string const & context = {} )
 {
   if( errorCode != 0 )
   {
     HYPREDRV_ErrorCodeDescribe( errorCode );
-    GEOS_ERROR( GEOS_FMT( "Error in call to {}", call ) );
+    std::fflush( stderr );
+    HYPREDRV_ErrorCodeClear();
+
+    std::string const contextMessage = context.empty()
+                                             ? std::string{}
+                                             : GEOS_FMT( "\nContext: {}", context );
+    GEOS_ERROR( GEOS_FMT( "Error in call to {} (HYPREDRV error code: 0x{:08x}, decimal: {}).{}\n"
+                          "HypreDrive error details are written to stderr; use '2>&1' to merge them into stdout. "
+                          "For additional library tracing, set HYPREDRV_LOG_LEVEL=3 and "
+                          "HYPREDRV_LOG_STREAM=stdout.",
+                          call,
+                          errorCode,
+                          errorCode,
+                          contextMessage ) );
   }
 }
 
@@ -401,6 +417,7 @@ bool strategyUsesCompositionalSemanticLabels( LinearSolverParameters::MGR::Strat
     case StrategyType::hybridSinglePhasePoromechanics:
     case StrategyType::singlePhasePoromechanicsEmbeddedFractures:
     case StrategyType::singlePhasePoromechanicsConformingFractures:
+    case StrategyType::singlePhasePoromechanicsConformingFracturesALM:
     case StrategyType::singlePhasePoromechanicsReservoirFVM:
     case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM:
     case StrategyType::hydrofracture:
@@ -869,7 +886,8 @@ enum class AMGFlavor
   pressure,
   pressureTemperature,
   displacementFiltered,
-  displacement
+  displacement,
+  almDisplacement
 };
 
 struct LevelAMGBlock
@@ -906,6 +924,7 @@ struct MGRSpecialization
   stdVector< LevelAMGBlock > fRelaxAMGLevels;
   HYPRE_Int pmax = 0;
   HYPRE_Int coarseMinCoarseSize = -1;
+  char const * cycle = nullptr;
 };
 
 MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType const strategy )
@@ -952,10 +971,12 @@ MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType c
     }
     case StrategyType::singlePhasePoromechanicsEmbeddedFractures:
     case StrategyType::singlePhasePoromechanicsConformingFractures:
+    case StrategyType::singlePhasePoromechanicsConformingFracturesALM:
     {
       MGRSpecialization specialization;
       specialization.coarseFlavor = AMGFlavor::pressure;
-      specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::displacement } };
+      specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::almDisplacement } };
+      specialization.cycle = "v(1,0)";
       return specialization;
     }
     case StrategyType::invalid:
@@ -1002,7 +1023,8 @@ void appendAMGHeader( std::ostringstream & stream,
 void appendDisplacementAMG( std::ostringstream & stream,
                             integer const indentLevel,
                             integer const separateComponents,
-                            bool const filterFunctions )
+                            bool const filterFunctions,
+                            bool const useALMSmoother = false )
 {
   appendAMGHeader( stream, indentLevel );
   appendLine( stream, indentLevel + 1, "coarsening:" );
@@ -1019,7 +1041,18 @@ void appendDisplacementAMG( std::ostringstream & stream,
   appendLine( stream, indentLevel + 2, "num_sweeps: 1" );
 #else
   appendLine( stream, indentLevel + 1, "relaxation:" );
-  appendLine( stream, indentLevel + 2, "order: 1" );
+  if( useALMSmoother )
+  {
+    appendLine( stream, indentLevel + 2, "down_type: l1sym-hgs" );
+    appendLine( stream, indentLevel + 2, "up_type: l1sym-hgs" );
+    appendLine( stream, indentLevel + 2, "coarse_type: ge" );
+    appendLine( stream, indentLevel + 2, "num_sweeps: 1" );
+    appendLine( stream, indentLevel + 2, "order: 0" );
+  }
+  else
+  {
+    appendLine( stream, indentLevel + 2, "order: 1" );
+  }
 #endif
 }
 
@@ -1115,6 +1148,11 @@ void appendAMGByFlavor( std::ostringstream & stream,
       appendDisplacementAMG( stream, indentLevel, 0, false );
       break;
     }
+    case AMGFlavor::almDisplacement:
+    {
+      appendDisplacementAMG( stream, indentLevel, mgrParams.separateComponents, true, true );
+      break;
+    }
   }
 }
 
@@ -1164,6 +1202,10 @@ bool buildStrategyYaml( LinearSolverParameters const & params,
   appendLine( stream, 2, "tolerance: 0.0" );
   appendLine( stream, 2, "max_iter: 1" );
   appendLine( stream, 2, GEOS_FMT( "print_level: {}", getMGRPrintLevel( params.logLevel ) ) );
+  if( specialization.cycle != nullptr )
+  {
+    appendLine( stream, 2, GEOS_FMT( "cycle: {}", specialization.cycle ) );
+  }
   appendLine( stream, 2, GEOS_FMT( "non_c_to_f: {}", 1 ) );
   appendLine( stream, 2, GEOS_FMT( "nonglk_max_elmts: {}", 1 ) );
   appendLine( stream, 2, GEOS_FMT( "pmax: {}", specialization.pmax ) );
@@ -1269,6 +1311,8 @@ bool buildMGRPreconditionerYaml( LinearSolverParameters const & params,
       return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsEmbeddedFractures >( params, labelNames, numComponentsPerField, preconditionerYaml );
     case StrategyType::singlePhasePoromechanicsConformingFractures:
       return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFractures >( params, labelNames, numComponentsPerField, preconditionerYaml );
+    case StrategyType::singlePhasePoromechanicsConformingFracturesALM:
+      return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFracturesALM >( params, labelNames, numComponentsPerField, preconditionerYaml );
     case StrategyType::singlePhasePoromechanicsReservoirFVM:
       return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsReservoirFVM >( params, labelNames, numComponentsPerField, preconditionerYaml );
     case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM:
@@ -1580,6 +1624,11 @@ void HypredriveSolver::setExecutionContext( LinearSolverExecutionContext const &
   m_hasExecutionContext = true;
 }
 
+void HypredriveSolver::setNearNullKernel( arrayView1d< HypreVector const > const & nearNullKernel )
+{
+  m_nearNullKernel = nearNullKernel;
+}
+
 void HypredriveSolver::setup( HypreMatrix const & mat )
 {
   Base::setup( mat );
@@ -1603,8 +1652,13 @@ void HypredriveSolver::createHypredrive( HypreMatrix const & mat,
   checkHypredriveCall( HYPREDRV_SetLibraryMode( m_hypredrive ), "HYPREDRV_SetLibraryMode" );
 
   char * argv[] = { const_cast< char * >( parseTarget.argument.c_str() ) };
+  std::string const parseContext =
+    parseTarget.source == hypre::hypredrive::InputSource::authoritativeFile
+    ? GEOS_FMT( "authoritative YAML file '{}'", parseTarget.argument )
+    : "YAML generated by GEOS";
   checkHypredriveCall( HYPREDRV_InputArgsParse( 1, argv, m_hypredrive ),
-                       "HYPREDRV_InputArgsParse" );
+                       "HYPREDRV_InputArgsParse",
+                       parseContext );
   if( parseTarget.source == hypre::hypredrive::InputSource::generatedFallback &&
       m_hasExecutionContext &&
       !m_executionContext.solverName.empty() )
@@ -1650,6 +1704,32 @@ void HypredriveSolver::refreshBoundObjects( HypreMatrix const & mat,
                                                          pointMarkers.data() ),
                          "HYPREDRV_LinearSystemSetDofmap" );
   }
+
+  if( !m_nearNullKernel.empty() )
+  {
+    localIndex const numEntries = mat.numLocalRows();
+    localIndex const numModes = m_nearNullKernel.size();
+    std::vector< HYPRE_Complex > values( numEntries * numModes );
+
+    for( localIndex mode = 0; mode < numModes; ++mode )
+    {
+      GEOS_ERROR_IF( m_nearNullKernel[mode].localSize() != numEntries,
+                     "HypreDrive near-null-space vector size does not match the matrix local size" );
+      GEOS_LAI_CHECK_ERROR(
+        HYPRE_IJVectorGetValues( m_nearNullKernel[mode].unwrappedIJ(),
+                                 LvArray::integerConversion< HYPRE_Int >( numEntries ),
+                                 nullptr,
+                                 values.data() + mode * numEntries ) );
+    }
+
+    checkHypredriveCall(
+      HYPREDRV_LinearSystemSetNearNullSpace(
+        m_hypredrive,
+        LvArray::integerConversion< int >( numEntries ),
+        LvArray::integerConversion< int >( numModes ),
+        values.data() ),
+      "HYPREDRV_LinearSystemSetNearNullSpace" );
+  }
 }
 
 void HypredriveSolver::setupLegacy( HypreMatrix const & mat )

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

@@ -98,6 +98,12 @@ class HypredriveSolver final : public LinearSolverBase< HypreInterface >
    */
   void setExecutionContext( LinearSolverExecutionContext const & context ) override;
 
+  /**
+   * @brief Set near-null-space modes to pass to HypreDrive during setup.
+   * @param nearNullKernel Full-system distributed near-null-space vectors.
+   */
+  void setNearNullKernel( arrayView1d< HypreVector const > const & nearNullKernel ) override;
+
   /**
    * @brief Build or refresh the solver/preconditioner for a matrix.
    * @param mat Matrix that defines the system structure and coefficients.
@@ -164,6 +170,7 @@ class HypredriveSolver final : public LinearSolverBase< HypreInterface >
   bool m_hasExecutionContext = false;
   bool m_timestepScopeActive = false;
   bool m_newtonScopeActive = false;
+  arrayView1d< HypreVector const > m_nearNullKernel;
   size_t m_hypredriveGeneration = 0;
   HYPREDRV_t m_hypredrive{};
   mutable HypreVector m_dummyRhs;