Merge branch 'main' into allowSolversToModify_621

CITATION.cff

@@ -1,7 +1,7 @@
 cff-version: 1.2.0
 message: If you use this software, please cite it as below.
 title: Model Independent Chemistry Model (MICM)
-version: v3.5.0
+version: v3.6.0
 doi: "10.5281/zenodo.10472189"
 authors:
   - family-names: Dawson

CMakeLists.txt

@@ -3,7 +3,7 @@
 cmake_minimum_required(VERSION 3.21)
 
 # project and version must be on the same line so that the docs can extract it
-project(micm VERSION 3.5.0 LANGUAGES CXX)
+project(micm VERSION 3.6.0 LANGUAGES CXX)
 
 if(NOT CMAKE_BUILD_TYPE)
   set(CMAKE_BUILD_TYPE "Release" CACHE STRING

docs/source/_static/switcher.json

@@ -1,7 +1,7 @@
 [
   {
-      "name": "v3.5.0 (stable)",
-      "version": "v3.5.0 (stable)",
+      "name": "v3.6.0 (stable)",
+      "version": "v3.6.0 (stable)",
       "url": "https://ncar.github.io/micm"
   },
   {
@@ -29,6 +29,11 @@
     "version": "3.5.0",
     "url": "https://ncar.github.io/micm/versions/3.5.0"
   },
+  {
+    "name": "v3.6.0",
+    "version": "3.6.0",
+    "url": "https://ncar.github.io/micm/versions/3.6.0"
+  },
   {
     "name": "dev",
     "version": "dev",

include/micm/cuda/solver/cuda_lu_decomposition.cuh

@@ -14,8 +14,7 @@ namespace micm
         const CudaMatrixParam& A_param,
         CudaMatrixParam& L_param,
         CudaMatrixParam& U_param,
-        const LuDecomposeParam& devstruct,
-        bool& is_singular);
+        const LuDecomposeParam& devstruct);
 
     /// This is the function that will copy the constant data
     ///   members of class "CudaLuDecomposition" to the device;

include/micm/cuda/solver/cuda_lu_decomposition.hpp

@@ -85,42 +85,21 @@ namespace micm
     /// @param A is the sparse matrix to decompose
     /// @param L is the lower triangular matrix created by decomposition
     /// @param U is the upper triangular matrix created by decomposition
-    /// @param is_singular Flag that is set to true if A is singular; false otherwise
-    template<class SparseMatrixPolicy>
-    requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void Decompose(
-        const SparseMatrixPolicy& A,
-        SparseMatrixPolicy& L,
-        SparseMatrixPolicy& U,
-        bool& is_singular) const;
-
     template<class SparseMatrixPolicy>
     requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void Decompose(
         const SparseMatrixPolicy& A,
         SparseMatrixPolicy& L,
         SparseMatrixPolicy& U) const;
   };
 
-  template<class SparseMatrixPolicy>
-  requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void CudaLuDecomposition::Decompose(
-      const SparseMatrixPolicy& A,
-      SparseMatrixPolicy& L,
-      SparseMatrixPolicy& U,
-      bool& is_singular) const
-  {
-    auto L_param = L.AsDeviceParam();  // we need to update lower matrix so it can't be constant and must be an lvalue
-    auto U_param = U.AsDeviceParam();  // we need to update upper matrix so it can't be constant and must be an lvalue
-    micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_, is_singular);
-  }
-
   template<class SparseMatrixPolicy>
   requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void CudaLuDecomposition::Decompose(
       const SparseMatrixPolicy& A,
       SparseMatrixPolicy& L,
       SparseMatrixPolicy& U) const
   {
-    bool is_singular = false;
     auto L_param = L.AsDeviceParam();  // we need to update lower matrix so it can't be constant and must be an lvalue
     auto U_param = U.AsDeviceParam();  // we need to update upper matrix so it can't be constant and must be an lvalue
-    micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_, is_singular);
+    micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_);
   }
 }  // end of namespace micm

include/micm/cuda/util/cuda_param.hpp

@@ -33,7 +33,6 @@ struct ProcessSetParam
 /// This struct could be allocated on the host or device;
 struct LuDecomposeParam
 {
-  bool* is_singular = nullptr;
   std::pair<std::size_t, std::size_t>* niLU_ = nullptr;
   char* do_aik_ = nullptr;
   std::size_t* aik_ = nullptr;

include/micm/jit/solver/jit_linear_solver.hpp

@@ -42,12 +42,7 @@ namespace micm
 
     /// @brief Decompose the matrix into upper and lower triangular matrices and general JIT functions
     /// @param matrix Matrix that will be factored into lower and upper triangular matrices
-    /// @param is_singular Flag that will be set to true if matrix is singular; false otherwise
-    void Factor(
-        SparseMatrixPolicy& matrix,
-        SparseMatrixPolicy& lower_matrix,
-        SparseMatrixPolicy& upper_matrix,
-        bool& is_singular) const;
+    void Factor(SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix) const;
 
     /// @brief Solve for x in Ax = b
     template<class MatrixPolicy>

include/micm/jit/solver/jit_linear_solver.inl

@@ -60,10 +60,9 @@ namespace micm
   inline void JitLinearSolver<L, SparseMatrixPolicy, LuDecompositionPolicy>::Factor(
       SparseMatrixPolicy &matrix,
       SparseMatrixPolicy &lower_matrix,
-      SparseMatrixPolicy &upper_matrix,
-      bool &is_singular) const
+      SparseMatrixPolicy &upper_matrix) const
   {
-    LinearSolver<SparseMatrixPolicy, LuDecompositionPolicy>::Factor(matrix, lower_matrix, upper_matrix, is_singular);
+    LinearSolver<SparseMatrixPolicy, LuDecompositionPolicy>::Factor(matrix, lower_matrix, upper_matrix);
   }
 
   template<std::size_t L, class SparseMatrixPolicy, class LuDecompositionPolicy>

include/micm/jit/solver/jit_lu_decomposition.hpp

@@ -39,10 +39,8 @@ namespace micm
     /// @param A Sparse matrix that will be decomposed
     /// @param lower The lower triangular matrix created by decomposition
     /// @param upper The upper triangular matrix created by decomposition
-    /// @param is_singular Flag that will be set to true if A is singular; false otherwise
     template<class SparseMatrixPolicy>
-    void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper, bool &is_singular)
-        const;
+    void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper) const;
 
    private:
     /// @brief Generates a function to perform the LU decomposition for a specific matrix sparsity structure

include/micm/jit/solver/jit_lu_decomposition.inl

@@ -89,6 +89,16 @@ namespace micm
           func.SetArrayElement(func.arguments_[2], U_ptr_index, JitType::Double, A_val);
           func.EndLoop(loop);
         }
+        else
+        {
+          auto loop = func.StartLoop("Uik_eq_zero_loop", 0, L);
+          llvm::Value *zero_val = llvm::ConstantFP::get(*(func.context_), llvm::APFloat(0.0));
+          llvm::Value *iUf = llvm::ConstantInt::get(*(func.context_), llvm::APInt(64, uik_nkj->first));
+          llvm::Value *U_ptr_index[1];
+          U_ptr_index[0] = func.builder_->CreateNSWAdd(loop.index_, iUf);
+          func.SetArrayElement(func.arguments_[2], U_ptr_index, JitType::Double, zero_val);
+          func.EndLoop(loop);
+        }
         for (std::size_t ikj = 0; ikj < uik_nkj->second; ++ikj)
         {
           auto loop = func.StartLoop("Uik_seq_Lij_Ujk_loop", 0, L);
@@ -137,6 +147,16 @@ namespace micm
           func.SetArrayElement(func.arguments_[1], L_ptr_index, JitType::Double, A_val);
           func.EndLoop(loop);
         }
+        else
+        {
+          auto loop = func.StartLoop("Lki_eq_zero_loop", 0, L);
+          llvm::Value *zero_val = llvm::ConstantFP::get(*(func.context_), llvm::APFloat(0.0));
+          llvm::Value *iLf = llvm::ConstantInt::get(*(func.context_), llvm::APInt(64, lki_nkj->first));
+          llvm::Value *L_ptr_index[1];
+          L_ptr_index[0] = func.builder_->CreateNSWAdd(loop.index_, iLf);
+          func.SetArrayElement(func.arguments_[1], L_ptr_index, JitType::Double, zero_val);
+          func.EndLoop(loop);
+        }
         for (std::size_t ikj = 0; ikj < lki_nkj->second; ++ikj)
         {
           auto loop = func.StartLoop("Lki_seq_Lkj_Uji_loop", 0, L);
@@ -186,25 +206,13 @@ namespace micm
 
   template<std::size_t L>
   template<class SparseMatrixPolicy>
-  void JitLuDecomposition<L>::Decompose(
-      const SparseMatrixPolicy &A,
-      SparseMatrixPolicy &lower,
-      SparseMatrixPolicy &upper,
-      bool &is_singular) const
+  void JitLuDecomposition<L>::Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper)
+      const
   {
-    is_singular = false;
     decompose_function_(A.AsVector().data(), lower.AsVector().data(), upper.AsVector().data());
     for (size_t block = 0; block < A.NumberOfBlocks(); ++block)
     {
       auto diagonals = upper.DiagonalIndices(block);
-      for (const auto &diagonal : diagonals)
-      {
-        if (upper.AsVector()[diagonal] == 0)
-        {
-          is_singular = true;
-          return;
-        }
-      }
     }
   }
 }  // namespace micm

include/micm/solver/backward_euler.inl

@@ -66,8 +66,6 @@ namespace micm
     std::size_t n_successful_integrations = 0;
     std::size_t n_convergence_failures = 0;
 
-    bool singular = false;
-
     auto derived_class_temporary_variables =
         static_cast<BackwardEulerTemporaryVariables<MatrixPolicy>*>(state.temporary_variables_.get());
     auto& Yn = derived_class_temporary_variables->Yn_;
@@ -112,7 +110,7 @@ namespace micm
         // (y_{n+1} - y_n) / H = f(t_{n+1}, y_{n+1})
 
         // try to find the root by factoring and solving the linear system
-        linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular);
+        linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_);
         result.stats_.decompositions_++;
 
         // forcing_blk in camchem

include/micm/solver/linear_solver.hpp

@@ -78,12 +78,7 @@ namespace micm
 
     /// @brief Decompose the matrix into upper and lower triangular matrices
     /// @param matrix Matrix to decompose into lower and upper triangular matrices
-    /// @param is_singular Flag that is set to true if matrix is singular; false otherwise
-    void Factor(
-        const SparseMatrixPolicy& matrix,
-        SparseMatrixPolicy& lower_matrix,
-        SparseMatrixPolicy& upper_matrix,
-        bool& is_singular) const;
+    void Factor(const SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix) const;
 
     /// @brief Solve for x in Ax = b. x should be a copy of b and after Solve finishes x will contain the result
     template<class MatrixPolicy>

include/micm/solver/linear_solver.inl

@@ -111,12 +111,11 @@ namespace micm
   inline void LinearSolver<SparseMatrixPolicy, LuDecompositionPolicy>::Factor(
       const SparseMatrixPolicy& matrix,
       SparseMatrixPolicy& lower_matrix,
-      SparseMatrixPolicy& upper_matrix,
-      bool& is_singular) const
+      SparseMatrixPolicy& upper_matrix) const
   {
     MICM_PROFILE_FUNCTION();
 
-    lu_decomp_.template Decompose<SparseMatrixPolicy>(matrix, lower_matrix, upper_matrix, is_singular);
+    lu_decomp_.template Decompose<SparseMatrixPolicy>(matrix, lower_matrix, upper_matrix);
   }
 
   template<class SparseMatrixPolicy, class LuDecompositionPolicy>

include/micm/solver/lu_decomposition.hpp

@@ -43,6 +43,13 @@ namespace micm
   /// For the sparse matrix algorithm, the indices of non-zero terms are stored in
   /// several arrays during construction. These arrays are iterated through during
   /// calls to Decompose to do the actual decomposition.
+  /// Our LU Decomposition only assigns the values of the jacobian to the LU matrices
+  /// when the *jacobian* is nonzero. However, the sparsity pattern of the jacobian doesn't
+  /// necessarily match that of the LU matrices. There can be more nonzero elements in the LU matrices
+  /// than in the jacobian. When this happens, we still need to assign the value of the jacobian matrix
+  /// to the LU matrix. This value is implicitly zero when the sparsity pattern differs. The Fill values
+  /// here do this implicit assignment
+  /// More detail in this issue: https://github.com/NCAR/micm/issues/625
   class LuDecomposition
   {
    protected:
@@ -115,19 +122,16 @@ namespace micm
     /// @param A Sparse matrix to decompose
     /// @param L The lower triangular matrix created by decomposition
     /// @param U The upper triangular matrix created by decomposition
-    /// @param is_singular Flag that is set to true if A is singular; false otherwise
     template<class SparseMatrixPolicy>
     requires(!VectorizableSparse<SparseMatrixPolicy>) void Decompose(
         const SparseMatrixPolicy& A,
         SparseMatrixPolicy& L,
-        SparseMatrixPolicy& U,
-        bool& is_singular) const;
+        SparseMatrixPolicy& U) const;
     template<class SparseMatrixPolicy>
     requires(VectorizableSparse<SparseMatrixPolicy>) void Decompose(
         const SparseMatrixPolicy& A,
         SparseMatrixPolicy& L,
-        SparseMatrixPolicy& U,
-        bool& is_singular) const;
+        SparseMatrixPolicy& U) const;
 
    private:
     /// @brief Initialize arrays for the LU decomposition