Hypre solver eliminate boundaries (Redux)

CMakeLists.txt

@@ -371,6 +371,7 @@ set(BOUT_SOURCES
     ./src/sys/generator_context.cxx
     ./include/bout/hyprelib.hxx
     ./src/sys/hyprelib.cxx
+    ./src/sys/hypre_interface.cxx
     ./src/sys/msg_stack.cxx
     ./src/sys/options.cxx
     ./src/sys/options/optionparser.hxx

include/bout/hypre_interface.hxx

@@ -5,8 +5,11 @@
 
 #if BOUT_HAS_HYPRE
 
+#include "bout/assert.hxx"
 #include "bout/bout_enum_class.hxx"
+#include "bout/bout_types.hxx"
 #include "bout/boutcomm.hxx"
+#include "bout/boutexception.hxx"
 #include "bout/caliper_wrapper.hxx"
 #include "bout/field.hxx"
 #include "bout/globalindexer.hxx"
@@ -58,6 +61,110 @@ int checkHypreError(int error) {
 // TODO: set sizes
 // TODO: set contiguous blocks at once
 
+/// Wrapper around HYPRE_Complex, that calls HypreFree when destroyed.
+struct HypreComplexArray {
+  HYPRE_Complex* data;
+
+  HypreComplexArray(int n) { HypreMalloc(data, sizeof(HYPRE_Complex) * n); }
+
+  ~HypreComplexArray() { HypreFree(data); }
+};
+
+/// Shared pointter to a HypreComplexArray. When the last copy is destroyed
+/// the HYPRE_Complex array inside will be free'd.
+using BCValuesPtr = std::shared_ptr<HypreComplexArray>;
+
+/*!
+ * This function modifies the input for the HYPRE_IJMatrixSetValues() routine to
+ * eliminate the boundary condition equations (see below for details on how the
+ * equations are adjusted).  It modifies the arrays ncols, rows, cols, and
+ * values.  It also returns a row_indexes array.  This can then be passed to the
+ * HYPRE_IJMatrixSetValues2() routine to set up the matrix in hypre.
+ *
+ * The arguments nb and bi_array indicate the boundary equations.  The routine
+ * returns info needed to adjust the right-hand-side and solution vector through
+ * the functions AdjustRightHandSideEquations and AdjustSolutionEquations.
+ *
+ * NOTE: It may make sense from an organizational standpoint to collect many of
+ * these arguments in a structure of some sort.
+ *
+ * Notation, assumptions, and other details:
+ *
+ * - Boundary equation i is assumed to have two coefficients
+ *
+ *      b_ii * u_i + b_ij * u_j = rhs_i
+ *
+ * - We also assume that each boundary equation has only one interior equation k
+ *   coupled to it (such that k = j) with coupling coefficient a_ki
+ *
+ *      a_ki * u_i + a_kj * u_j + ... = rhs_k
+ *
+ * - Each equation k is adjusted as follows:
+ *
+ *      a_kj = a_kj - a_ki * b_ij / b_ii
+ *      a_ki = 0
+ *
+ * - Boundary equations are adjusted to be identity equations in the matrix, but
+ *   the boundary coefficients (b_ii, b_ij) are returned for use later
+ *
+ * - Right-hand-side equations are adjusted in AdjustRightHandSideEquations() as
+ *   follows: rhs_k = rhs_k - a_ki * rhs_i / b_ii
+ *
+ * - Solution unknowns are adjusted at boundaries in AdjustSolutionEquations as
+ *   follows: u_i = (rhs_i - b_ij * u_j) / b_ii
+ *
+ * - Naming conventions: Arrays starting with 'b' are boundary equation arrays
+ *   indexed by 'bnum', and arrays starting with 'a' are non-boundary arrays
+ *   (interior matrix equations) indexed by 'anum'.  When 'num' is prefixed with
+ *   a row or column number 'i', 'j', or 'k', the array holds the corresponding
+ *   local data index for that row or column (e.g., an index into the local
+ *   solution vector).  Matrix coefficients are named as above, e.g., 'bij' is
+ *   the coefficient for b_ij.
+ *
+ *   NOTE: Implementation in src/sys/hypre_interface.cxx
+ */
+struct BCMatrixEquations {
+  HYPRE_Int nb;
+  HYPRE_Int* binum_array;
+  HYPRE_Int* bjnum_array;
+  HYPRE_Complex* bii_array;
+  HYPRE_Complex* bij_array;
+  HYPRE_Int na;
+  HYPRE_Int* aknum_array;
+  HYPRE_Complex* aki_array;
+
+  BCMatrixEquations() = delete;
+
+  BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_BigInt* rows,
+                    HYPRE_Int** row_indexes_ptr, HYPRE_BigInt* cols,
+                    HYPRE_Complex* values,
+                    HYPRE_Int nb,         // number of boundary equations
+                    HYPRE_Int* bi_array); // row i for each boundary equation
+
+  ~BCMatrixEquations() {
+    // Free arrays
+    HypreFree(binum_array);
+    HypreFree(bjnum_array);
+    HypreFree(bii_array);
+    HypreFree(bij_array);
+    HypreFree(aknum_array);
+    HypreFree(aki_array);
+  }
+
+  /// Applies in-place modification of the rhs array.
+  ///
+  /// Returns an array of boundary values that can be used to apply
+  /// boundary conditions to a solution vector.
+  BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex* rhs);
+
+  /// Apply boundary conditions to the solution.
+  /// Uses the BCValuesPtr returned from adjustBCRightHandSideEquations()
+  void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex* solution);
+};
+
+/// A shared pointer to a BCMatrixEquations object
+using BCMatrixPtr = std::shared_ptr<BCMatrixEquations>;
+
 template <class T>
 class HypreVector {
   MPI_Comm comm;
@@ -174,6 +281,14 @@ public:
     HypreMalloc(V, vsize * sizeof(HYPRE_Complex));
   }
 
+  // Data for eliminating boundary equation
+  bool elimBErhs = false;
+  bool elimBEsol = false;
+  BCMatrixPtr bcmatrix;
+  BCValuesPtr bcvalues; /// Stores rhs values of BC rows
+
+  void syncElimBErhs(HypreVector<T>& rhs) { bcvalues = rhs.bcvalues; }
+
   void assemble() {
     CALI_CXX_MARK_FUNCTION;
     writeCacheToHypre();
@@ -183,11 +298,17 @@ public:
   }
 
   void writeCacheToHypre() {
+    if (elimBErhs) {
+      bcvalues = bcmatrix->adjustBCRightHandSideEquations(V);
+    }
     checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V));
   }
 
   void readCacheFromHypre() {
     checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V));
+    if (elimBEsol) {
+      bcmatrix->adjustBCSolutionEquations(bcvalues, V);
+    }
   }
 
   T toField() {
@@ -210,6 +331,7 @@ public:
     return result;
   }
 
+  // Loads values into I and V arrays but does not assemble
   void importValuesFromField(const T& f) {
     CALI_CXX_MARK_FUNCTION;
 
@@ -224,8 +346,6 @@ public:
     }
 
     ASSERT2(vec_i == vsize);
-    // writeCacheToHypre(); // redundant assemble already performs writeCacheToHypre
-    assemble();
     have_indices = true;
   }
 
@@ -265,19 +385,19 @@ public:
     }
 
     Element& operator=(const Element& other) {
-      ASSERT3(finite(static_cast<BoutReal>(other)));
+      ASSERT3(std::isfinite(static_cast<BoutReal>(other)));
       return *this = static_cast<BoutReal>(other);
     }
     Element& operator=(BoutReal value_) {
-      ASSERT3(finite(value_));
+      ASSERT3(std::isfinite(value_));
       value = value_;
       vector->V[vec_i] = value_;
       return *this;
     }
     Element& operator+=(BoutReal value_) {
-      ASSERT3(finite(value_));
+      ASSERT3(std::isfinite(value_));
       value += value_;
-      ASSERT3(finite(value));
+      ASSERT3(std::isfinite(value));
       vector->V[vec_i] += value_;
       return *this;
     }
@@ -432,7 +552,7 @@ public:
           weights(weights_) {
 #if CHECK > 2
       for (const auto val : weights) {
-        ASSERT3(finite(val));
+        ASSERT3(std::isfinite(val));
       }
 #endif
       ASSERT2(positions.size() == weights.size());
@@ -445,19 +565,19 @@ public:
     }
     Element& operator=(const Element& other) {
 
-      ASSERT3(finite(static_cast<BoutReal>(other)));
+      ASSERT3(std::isfinite(static_cast<BoutReal>(other)));
       return *this = static_cast<BoutReal>(other);
     }
     Element& operator=(BoutReal value_) {
 
-      ASSERT3(finite(value_));
+      ASSERT3(std::isfinite(value_));
       value = value_;
       setValues(value);
       return *this;
     }
     Element& operator+=(BoutReal value_) {
 
-      ASSERT3(finite(value_));
+      ASSERT3(std::isfinite(value_));
       auto column_position = std::find(cbegin(positions), cend(positions), column);
       if (column_position != cend(positions)) {
         const auto i = std::distance(cbegin(positions), column_position);
@@ -667,6 +787,20 @@ public:
     return Element(*this, global_row, global_column, positions, weights);
   }
 
+  // Data for eliminating boundary equations
+  bool elimBE = false;
+  BCMatrixPtr bcmatrix; // Shared pointer
+
+  void setElimBE() { elimBE = true; }
+
+  void setElimBEVectors(HypreVector<T>& sol, HypreVector<T>& rhs) {
+    sol.elimBEsol = elimBE;
+    sol.bcmatrix = bcmatrix;
+
+    rhs.elimBErhs = elimBE;
+    rhs.bcmatrix = bcmatrix;
+  }
+
   void assemble() {
     CALI_CXX_MARK_FUNCTION;
 
@@ -695,8 +829,32 @@ public:
         entry++;
       }
     }
-    checkHypreError(
-        HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals));
+
+    // Eliminate boundary condition equations in hypre SetValues input arguments
+    if (elimBE) {
+      HYPRE_Int* bi_array;
+      HYPRE_Int* row_indexes;
+      // There must be an easier way to get nb
+      int nb = 0;
+      BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { nb++; }
+      HypreMalloc(bi_array, nb * sizeof(HYPRE_Int));
+      nb = 0;
+      BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) {
+        bi_array[nb] = index_converter->getGlobal(i);
+        nb++;
+      }
+
+      bcmatrix = std::make_shared<BCMatrixEquations>(
+          num_rows, num_cols, rawI, &row_indexes, cols, vals, nb, bi_array);
+      HypreFree(bi_array);
+
+      checkHypreError(HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI,
+                                               row_indexes, cols, vals));
+      HypreFree(row_indexes);
+    } else {
+      checkHypreError(
+          HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals));
+    }
     checkHypreError(HYPRE_IJMatrixAssemble(*hypre_matrix));
     checkHypreError(HYPRE_IJMatrixGetObject(*hypre_matrix,
                                             reinterpret_cast<void**>(&parallel_matrix)));
@@ -876,6 +1034,31 @@ public:
     setMaxIter(
         options["maxits"].doc("Maximum iterations for Hypre solver").withDefault(10000));
 
+    switch (solver_type) {
+    case HYPRE_SOLVER_TYPE::gmres: {
+      HYPRE_ParCSRGMRESSetKDim(solver,
+                               options["kdim"]
+                                   .doc("Set the maximum size of the Krylov space")
+                                   .withDefault(30));
+
+      if (options["skip_real_residual_check"]
+              .doc("Skip the evaluation and the check of the actual residual?")
+              .withDefault<bool>(false)) {
+        HYPRE_GMRESSetSkipRealResidualCheck(solver, 1);
+      }
+      break;
+    }
+    case HYPRE_SOLVER_TYPE::bicgstab: {
+      break;
+    }
+    case HYPRE_SOLVER_TYPE::pcg: {
+      break;
+    }
+    default: {
+      throw BoutException("Unsupported hypre_solver_type {}", toString(solver_type));
+    }
+    }
+
     HYPRE_BoomerAMGCreate(&precon);
     HYPRE_BoomerAMGSetOldDefault(precon);
 #if BOUT_HAS_CUDA

src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx

@@ -218,11 +218,15 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) {
 
   CALI_MARK_BEGIN("LaplaceHypre3d_solve:vectorAssemble");
 
+  operator3D.setElimBEVectors(solution, rhs);
+
   rhs.importValuesFromField(b);
   solution.importValuesFromField(x0);
   rhs.assemble();
   solution.assemble();
 
+  solution.syncElimBErhs(rhs);
+
   CALI_MARK_END("LaplaceHypre3d_solve:vectorAssemble");
 
   CALI_MARK_BEGIN("LaplaceHypre3d_solve:solve");
@@ -411,6 +415,7 @@ void LaplaceHypre3d::updateMatrix3D() {
     operator3D.ydown(ydown)(l, l.ym().zp()) += -C_d2f_dydz;
     operator3D.ydown(ydown)(l, l.ym().zm()) += C_d2f_dydz;
   }
+  operator3D.setElimBE();
   operator3D.assemble();
 
   if (print_matrix) {