572 check for singularity when the solver parameters flag is turned on

include/micm/solver/lu_decomposition.hpp

@@ -108,16 +108,6 @@ namespace micm
         const SparseMatrixPolicy& A,
         typename SparseMatrixPolicy::value_type initial_value);
 
-    /// @brief Perform an LU decomposition on a given A matrix
-    /// @param A Sparse matrix to decompose
-    /// @param L The lower triangular matrix created by decomposition
-    /// @param U The upper triangular matrix created by decomposition
-    template<class SparseMatrixPolicy>
-    requires(SparseMatrixConcept<SparseMatrixPolicy>) void Decompose(
-        const SparseMatrixPolicy& A,
-        SparseMatrixPolicy& L,
-        SparseMatrixPolicy& U) const;
-
     /// @brief Perform an LU decomposition on a given A matrix
     /// @param A Sparse matrix to decompose
     /// @param L The lower triangular matrix created by decomposition

include/micm/solver/lu_decomposition.inl

@@ -163,16 +163,6 @@ namespace micm
     return LU;
   }
 
-  template<class SparseMatrixPolicy>
-  requires(SparseMatrixConcept<SparseMatrixPolicy>) inline void LuDecomposition::Decompose(
-      const SparseMatrixPolicy& A,
-      SparseMatrixPolicy& L,
-      SparseMatrixPolicy& U) const
-  {
-    bool is_singular = false;
-    Decompose<SparseMatrixPolicy>(A, L, U, is_singular);
-  }
-
   template<class SparseMatrixPolicy>
   requires(!VectorizableSparse<SparseMatrixPolicy>) inline void LuDecomposition::Decompose(
       const SparseMatrixPolicy& A,
@@ -223,16 +213,19 @@ namespace micm
             L_vector[lki_nkj->first] -= L_vector[lkj_uji->first] * U_vector[lkj_uji->second];
             ++lkj_uji;
           }
+
           if (U_vector[*uii] == 0.0)
           {
             is_singular = true;
-            return;
           }
           L_vector[lki_nkj->first] /= U_vector[*uii];
           ++lki_nkj;
           ++uii;
         }
       }
+      // the singularity check inside the for loop won't detect a zero in the bottom right position
+      auto cell_U_bottom_right = std::next(U.AsVector().begin(), i_block * U.FlatBlockSize() + U.FlatBlockSize() - 1);
+      if (*cell_U_bottom_right == 0) is_singular = true;
     }
   }
 
@@ -316,15 +309,17 @@ namespace micm
             if (U_vector[uii_deref + i_cell] == 0.0)
             {
               is_singular = true;
-              return;
             }
             L_vector[lki_nkj_first + i_cell] /= U_vector[uii_deref + i_cell];
           }
           ++lki_nkj;
           ++uii;
         }
       }
+
+      auto cell_U_bottom_right = std::next(U.AsVector().begin(), i_group * U_GroupSizeOfFlatBlockSize + U_GroupSizeOfFlatBlockSize - 1);
+      if (*cell_U_bottom_right == 0) is_singular = true;
     }
   }
 
-}  // namespace micm
+}  // namespace micm

include/micm/solver/rosenbrock.inl

@@ -240,26 +240,21 @@ namespace micm
   {
     MICM_PROFILE_FUNCTION();
 
-    auto jacobian = state.jacobian_;
     uint64_t n_consecutive = 0;
     singular = false;
     while (true)
     {
+      auto jacobian = state.jacobian_;
       double alpha = 1 / (H * gamma);
       static_cast<Derived*>(this)->AlphaMinusJacobian(jacobian, alpha);
-      if (parameters_.check_singularity_)
-      {
-        linear_solver_.Factor(jacobian, state.lower_matrix_, state.upper_matrix_, singular);
-      }
-      else
-      {
-        singular = false;
-        linear_solver_.Factor(jacobian, state.lower_matrix_, state.upper_matrix_);
-      }
-      singular = false;  // TODO This should be evaluated in some way
+      linear_solver_.Factor(jacobian, state.lower_matrix_, state.upper_matrix_, singular);
       stats.decompositions_ += 1;
-      if (!singular)
+
+      // if we are checking for singularity and the matrix is not singular, we can break the loop
+      // if we are not checking for singularity, we always break the loop
+      if (!singular || !parameters_.check_singularity_)
         break;
+
       stats.singular_ += 1;
       if (++n_consecutive > 5)
         break;

include/micm/solver/state.inl

@@ -92,7 +92,7 @@ namespace micm
     jacobian_ = BuildJacobian<SparseMatrixPolicy>(
         parameters.nonzero_jacobian_elements_, parameters.number_of_grid_cells_, state_size_);
 
-    auto lu = LuDecomposition::GetLUMatrices<SparseMatrixPolicy>(jacobian_, 1.0e-30);
+    auto lu = LuDecomposition::GetLUMatrices<SparseMatrixPolicy>(jacobian_, 0);
     auto lower_matrix = std::move(lu.first);
     auto upper_matrix = std::move(lu.second);
     lower_matrix_ = lower_matrix;

test/integration/test_analytical_rosenbrock.cpp

@@ -220,7 +220,7 @@ TEST(AnalyticalExamples, SurfaceRxn)
   test_analytical_surface_rxn(rosenbrock_6stage_da, 1e-7);
 }
 
-TEST(AnalyticalExamples, HIRESConfig)
+TEST(AnalyticalExamples, HIRES)
 {
   auto rosenbrock_solver = [](auto params)
   {

test/unit/solver/test_rosenbrock.cpp

@@ -46,6 +46,70 @@ SolverBuilderPolicy getSolver(SolverBuilderPolicy builder)
       .SetReorderState(false);
 }
 
+template<class SolverBuilderPolicy>
+SolverBuilderPolicy getSingularSystemZeroInBottomRightOfU(SolverBuilderPolicy builder)
+{
+  // A -> B
+  // B -> A
+
+  auto a = micm::Species("a");
+  auto b = micm::Species("b");
+
+  micm::Phase gas_phase{ std::vector<micm::Species>{ a, b } };
+
+  micm::Process r1 = micm::Process::Create()
+                         .SetReactants({ a })
+                         .SetProducts({ Yields(b, 1) })
+                         .SetPhase(gas_phase)
+                         .SetRateConstant(micm::UserDefinedRateConstant({.label_ = "r1"}));
+
+  micm::Process r2 = micm::Process::Create()
+                          .SetReactants({ b })
+                          .SetProducts({ Yields(a, 1) })
+                          .SetPhase(gas_phase)
+                          .SetRateConstant(micm::UserDefinedRateConstant({.label_ = "r2"}));
+
+  return builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }))
+      .SetReactions(std::vector<micm::Process>{ r1, r2 })
+      .SetReorderState(false);
+}
+
+template<class SolverBuilderPolicy>
+SolverBuilderPolicy getSolverForSingularSystemOnDiagonal(SolverBuilderPolicy builder)
+{
+  // A -> B, k1
+  // B -> C, k2
+  // C -> A, k3
+
+  auto a = micm::Species("a");
+  auto b = micm::Species("b");
+  auto c = micm::Species("c");
+
+  micm::Phase gas_phase{ std::vector<micm::Species>{ a, b, c } };
+
+  micm::Process r1 = micm::Process::Create()
+                         .SetReactants({ a })
+                         .SetProducts({ Yields(b, 1) })
+                         .SetPhase(gas_phase)
+                         .SetRateConstant(micm::UserDefinedRateConstant({.label_ = "r1"}));
+
+  micm::Process r2 = micm::Process::Create()
+                          .SetReactants({ b })
+                          .SetProducts({ Yields(c, 1) })
+                          .SetPhase(gas_phase)
+                          .SetRateConstant(micm::UserDefinedRateConstant({.label_ = "r2"}));
+
+  micm::Process r3 = micm::Process::Create()
+                          .SetReactants({ c })
+                          .SetProducts({ Yields(a, 1) })
+                          .SetPhase(gas_phase)
+                          .SetRateConstant(micm::UserDefinedRateConstant({.label_ = "r3"}));
+
+  return builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }))
+      .SetReactions(std::vector<micm::Process>{ r1, r2, r3 })
+      .SetReorderState(false);
+}
+
 template<class SolverBuilderPolicy>
 void testAlphaMinusJacobian(SolverBuilderPolicy builder, std::size_t number_of_grid_cells)
 {
@@ -225,4 +289,103 @@ TEST(RosenbrockSolver, VectorNormalizedError)
   testNormalizedErrorDiff(VectorBuilder<4>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 3);
   testNormalizedErrorDiff(VectorBuilder<8>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 5);
   testNormalizedErrorDiff(VectorBuilder<10>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 3);
-}
+}
+
+TEST(RosenbrockSolver, SingularSystemZeroInBottomRightOfU)
+{
+  auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters();
+  params.check_singularity_ = true;
+  auto standard = StandardBuilder(params);
+  auto vector = VectorBuilder<3>(params);
+
+  auto standard_solver = getSingularSystemZeroInBottomRightOfU(standard).SetNumberOfGridCells(1).Build();
+  auto vector_solver = getSingularSystemZeroInBottomRightOfU(vector).SetNumberOfGridCells(4).Build();
+
+  auto standard_state = standard_solver.GetState();
+  auto vector_state = vector_solver.GetState();
+
+  double k1 = -2;
+  double k2 = 1.0;
+
+  standard_state.SetCustomRateParameter("r1", k1);
+  standard_state.SetCustomRateParameter("r2", k2);
+
+  vector_state.SetCustomRateParameter("r1", {k1, k1, k1, k1});
+  vector_state.SetCustomRateParameter("r2", {k2, k2, k2, k2});
+
+  standard_state.variables_[0] = { 1.0, 1.0 };
+  vector_state.variables_[0] = { 1.0, 1.0 };
+  vector_state.variables_[1] = { 1.0, 1.0 };
+  vector_state.variables_[2] = { 1.0, 1.0 };
+  vector_state.variables_[3] = { 1.0, 1.0 };
+
+  // to get a jacobian with an LU factorization that contains a zero on the diagonal
+  // of U, we need det(alpha * I - jacobian) = 0
+  // for the system above, that means we have to have alpha + k1 + k2 = 0
+  // in this case, one of the reaction rates will be negative but it's good enough to 
+  // test the singularity check
+  // alpha is 1 / (H * gamma), where H is the time step and gamma is the gamma value from 
+  // the rosenbrock paramters
+  // so H needs to be 1 / ( (-k1 - k2) * gamma)
+  // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1
+  double H = 1 / ( (-k1 - k2) * params.gamma_[0]);
+  standard_solver.solver_.parameters_.h_start_ = H;
+
+  standard_solver.CalculateRateConstants(standard_state);
+  vector_solver.CalculateRateConstants(vector_state);
+
+  auto standard_result = standard_solver.Solve(2*H, standard_state);
+  EXPECT_NE(standard_result.stats_.singular_, 0);
+
+  auto vector_result = vector_solver.Solve(2*H, vector_state);
+  EXPECT_NE(vector_result.stats_.singular_, 0);
+}
+
+TEST(RosenbrockSolver, SingularSystemZeroAlongDiagonalNotBottomRight)
+{
+  auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters();
+
+  double k1 = -1.0;
+  double k2 = -1.0;
+  double k3 = 1.0;
+
+  // to get a jacobian with an LU factorization that contains a zero on the diagonal
+  // of U, we need det(alpha * I - jacobian) = 0
+  // for the system above, that means we have to set alpha = -k1, or alpha=-k2, or alpha=k3
+  double H = 1 / ( -k1* params.gamma_[0]);
+
+  params.check_singularity_ = true;
+  params.h_start_ = H;
+
+  auto standard = StandardBuilder(params);
+  auto vector = VectorBuilder<3>(params);
+
+  auto standard_solver = getSolverForSingularSystemOnDiagonal(standard).SetNumberOfGridCells(1).Build();
+  auto vector_solver = getSolverForSingularSystemOnDiagonal(vector).SetNumberOfGridCells(4).Build();
+
+  auto standard_state = standard_solver.GetState();
+  auto vector_state = vector_solver.GetState();
+
+  standard_state.SetCustomRateParameter("r1", k1);
+  standard_state.SetCustomRateParameter("r2", k2);
+  standard_state.SetCustomRateParameter("r3", k3);
+
+  vector_state.SetCustomRateParameter("r1", {k1, k1, k1, k1});
+  vector_state.SetCustomRateParameter("r2", {k2, k2, k2, k2});
+  vector_state.SetCustomRateParameter("r3", {k3, k3, k3, k3});
+
+  standard_state.variables_[0] = { 1.0, 1.0, 1.0 };
+  vector_state.variables_[0] =   { 1.0, 1.0, 1.0 };
+  vector_state.variables_[1] =   { 1.0, 1.0, 1.0 };
+  vector_state.variables_[2] =   { 1.0, 1.0, 1.0 };
+  vector_state.variables_[3] =   { 1.0, 1.0, 1.0 };
+
+  standard_solver.CalculateRateConstants(standard_state);
+  vector_solver.CalculateRateConstants(vector_state);
+
+  auto standard_result = standard_solver.Solve(2*H, standard_state);
+  EXPECT_NE(standard_result.stats_.singular_, 0);
+
+  auto vector_result = vector_solver.Solve(2*H, vector_state);
+  EXPECT_NE(vector_result.stats_.singular_, 0);
+}