Merge pull request #1301 from LLNL/tupek/trust_region_petsc_errors

Fix a bunch of small errors.

src/serac/numerics/dense_petsc.hpp

@@ -16,7 +16,13 @@ struct DenseMat {
 
   /// @brief constructor
   /// @param a matrix
-  DenseMat(const DenseMat& a) { MatDuplicate(a.A, MAT_COPY_VALUES, &A); }
+  DenseMat(const DenseMat& a)
+  {
+    MatDuplicate(a.A, MAT_COPY_VALUES, &A);
+    MatCopy(a.A, A, SAME_NONZERO_PATTERN);
+    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
+  }
 
   /// @brief destructor
   ~DenseMat() { MatDestroy(&A); }
@@ -59,6 +65,26 @@ struct DenseMat {
     MatView(A, PETSC_VIEWER_STDOUT_SELF);
   }
 
+  /// @brief check for nans
+  bool hasNan() const
+  {
+    auto [rows, cols] = size();
+    for (int i = 0; i < rows; ++i) {
+      for (int j = 0; j < cols; ++j) {
+        double val = (*this)(i, j);
+        if (val != val) return true;
+      }
+    }
+    return false;
+  }
+
+  /// @brief  reassemble petsc dense matrix after values have been modified
+  void reassemble()
+  {
+    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
+  }
+
   /// petsc matrix
   Mat A;
 };
@@ -88,6 +114,9 @@ DenseMat sym(const DenseMat& a)
       b.setValue(j, i, val);
     }
   }
+
+  b.reassemble();
+
   return b;
 }
 

src/serac/numerics/equation_solver.cpp

@@ -382,12 +382,32 @@ class TrustRegion : public mfem::NewtonSolver {
       H_directions.emplace_back(H_left.get());
     }
 
-    std::tie(directions, H_directions) = removeDependentDirections(directions, H_directions);
+    try {
+      std::tie(directions, H_directions) = removeDependentDirections(directions, H_directions);
+    } catch (const std::exception& e) {
+      if (print_options.warnings) {
+        mfem::out << "remove dependent directions failed with " << e.what() << std::endl;
+      }
+      return;
+    }
 
     mfem::Vector b(g);
     b *= -1;
-    auto [sol, leftvecs, leftvals, energy_change] =
-        solveSubspaceProblem(directions, H_directions, b, delta, num_leftmost);
+
+    mfem::Vector                               sol;
+    std::vector<std::shared_ptr<mfem::Vector>> leftvecs;
+    std::vector<double>                        leftvals;
+    double                                     energy_change;
+
+    try {
+      std::tie(sol, leftvecs, leftvals, energy_change) =
+          solveSubspaceProblem(directions, H_directions, b, delta, num_leftmost);
+    } catch (const std::exception& e) {
+      if (print_options.warnings) {
+        mfem::out << "subspace solve failed with " << e.what() << std::endl;
+      }
+      return;
+    }
 
     left_mosts.clear();
     for (auto& lv : leftvecs) {
@@ -398,8 +418,9 @@ class TrustRegion : public mfem::NewtonSolver {
     double subspace_energy = computeEnergy(g, hess_vec_func, sol);
 
     if (print_options.iterations || print_options.warnings) {
+      double leftval = leftvals.size() ? leftvals[0] : 1.0;
       mfem::out << "Energy using subspace solver from: " << base_energy << ", to: " << subspace_energy << " / "
-                << energy_change << ".  Min eig: " << leftvals[0] << std::endl;
+                << energy_change << ".  Min eig: " << leftval << std::endl;
     }
 
     if (subspace_energy < base_energy) {
@@ -727,7 +748,6 @@ class TrustRegion : public mfem::NewtonSolver {
         if (use_with_option1 || use_with_option2 || use_with_option3) {
           if (!have_computed_Hvs) {
             have_computed_Hvs = true;
-
             hess_vec_func(trResults.z, trResults.H_z);
             hess_vec_func(trResults.d_old, trResults.H_d_old);
             hess_vec_func(trResults.cauchy_point, trResults.H_cauchy_point);
@@ -755,9 +775,8 @@ class TrustRegion : public mfem::NewtonSolver {
         double realObjective = std::numeric_limits<double>::max();
         double normPred      = std::numeric_limits<double>::max();
         try {
-          normPred    = computeResidual(x_pred, r_pred);
-          double obj1 = 0.5 * (Dot(r, trResults.d) + Dot(r_pred, trResults.d)) - roundOffTol;
-
+          normPred      = computeResidual(x_pred, r_pred);
+          double obj1   = 0.5 * (Dot(r, trResults.d) + Dot(r_pred, trResults.d)) - roundOffTol;
           realObjective = obj1;
         } catch (const std::exception&) {
           realObjective = std::numeric_limits<double>::max();

src/serac/numerics/trust_region_solver.cpp

@@ -185,7 +185,9 @@ auto qr(const std::vector<const mfem::Vector*>& states)
   Mat R;
   int num_cols = static_cast<int>(states.size());
   MatCreateSeqDense(PETSC_COMM_SELF, num_cols, num_cols, NULL, &R);
-  BVOrthogonalize(Q, R);
+  auto error = BVOrthogonalize(Q, R);
+
+  if (error) throw PetscException("BVOrthogonalize failed.");
 
   return std::make_pair(Q, DenseMat(R));
 }
@@ -383,8 +385,18 @@ std::tuple<mfem::Vector, std::vector<std::shared_ptr<mfem::Vector>>, std::vector
   DenseMat sAs1 = dot(states, Astates);
   DenseMat sAs  = sym(sAs1);
 
+  if (sAs.hasNan()) {
+    throw PetscException("States in subspace solve contain NaNs.");
+    return std::make_tuple(b, std::vector<std::shared_ptr<mfem::Vector>>{}, std::vector<double>{}, 0);
+  }
+
   auto [Q_parallel, R] = qr(states);
 
+  if (R.hasNan()) {
+    throw PetscException("R from qr returning with a NaN.");
+    return std::make_tuple(b, std::vector<std::shared_ptr<mfem::Vector>>{}, std::vector<double>{}, 0);
+  }
+
   auto [rows, cols] = R.size();
   SLIC_ERROR_IF(rows != cols, "R matrix is not square in subspace problem solve\n");
 
@@ -434,7 +446,6 @@ std::tuple<mfem::Vector, std::vector<std::shared_ptr<mfem::Vector>>, std::vector
   BVDestroy(&Q_parallel);
   VecDestroy(&b_parallel);
   VecDestroy(&x_parallel);
-
   return std::make_tuple(sol, leftmosts, leftvals, energy);
 }
 
@@ -443,6 +454,7 @@ std::tuple<mfem::Vector, std::vector<std::shared_ptr<mfem::Vector>>, std::vector
 std::pair<std::vector<const mfem::Vector*>, std::vector<const mfem::Vector*>> removeDependentDirections(
     std::vector<const mfem::Vector*> directions, std::vector<const mfem::Vector*> A_directions)
 {
+  SERAC_MARK_FUNCTION;
   std::vector<double> norms;
   size_t              num_dirs = directions.size();
 

src/serac/numerics/trust_region_solver.hpp

@@ -25,6 +25,19 @@
 
 namespace serac {
 
+class PetscException : public std::exception {
+public:
+  /// constructor
+  PetscException(const std::string& message) : msg(message) {}
+
+  /// what is message
+  const char* what() const noexcept override { return msg.c_str(); }
+
+private:
+  /// message string
+  std::string msg;
+};
+
 /// @brief computes the global size of mfem::Vector
 int globalSize(const mfem::Vector& parallel_v, const MPI_Comm& comm);