Fix substepping issue for Rosenbrock solver with inline LU decomposit…

…ion (#717)

* revise Rosenbrock sovler with substepping for inline LU; add integration tests for CSC and inline LU decomp

* bug fix for compilation

* add back the last_alpha value for non-inline LU algorithm

* minor bug fix

* minor style fix

* clean up the code further

* add suggested flag for Windows with MSVC compiler

* debug message

* try another fix

* move the location

* reduce parallel jobs for clang-cl-11

* add big object file compiler option for clang compiler

* add flag for clang-cl only

* remove win32 for clang-cl-11

* clean up cmake file

.github/workflows/windows.yml

@@ -77,7 +77,7 @@ jobs:
     continue-on-error: true 
     strategy:
       matrix:
-        architecture: [Win32, x64]
+        architecture: [x64]
 
     steps:
       - uses: actions/checkout@v4

CMakeLists.txt

@@ -46,6 +46,12 @@ if (${CMAKE_HOST_SYSTEM_NAME} MATCHES "Linux" AND "${CMAKE_CXX_COMPILER_ID}" STR
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++")
 endif()
 
+# on Windows with MSVC, add the /bigobj flag to allow for large object files
+if (${CMAKE_HOST_SYSTEM_NAME} MATCHES "Windows" AND "${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC")
+  # If the compiler is MSVC or Clang on Windows, use /bigobj
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /bigobj")
+endif()
+
 ################################################################################
 # Dependencies
 

include/micm/solver/rosenbrock.hpp

@@ -104,7 +104,7 @@ namespace micm
     /// @param number_densities The number densities
     /// @param stats The solver stats
     /// @param state The state
-    void LinearFactor(const double alpha, const auto& number_densities, SolverStats& stats, auto& state) const;
+    void LinearFactor(const double alpha, SolverStats& stats, auto& state) const;
 
     /// @brief Computes the scaled norm of the vector errors
     /// @param y the original vector

include/micm/solver/rosenbrock.inl

@@ -10,13 +10,14 @@ namespace micm
       const RosenbrockSolverParameters& parameters) const noexcept
   {
     MICM_PROFILE_FUNCTION();
-    using MatrixPolicy = decltype(state.variables_);
+    using DenseMatrixPolicy = decltype(state.variables_);
+    using SparseMatrixPolicy = decltype(state.jacobian_);
 
     SolverResult result{};
     result.state_ = SolverState::Running;
     auto& Y = state.variables_;  // Y will hold the new solution at the end of the solve
     auto derived_class_temporary_variables =
-        static_cast<RosenbrockTemporaryVariables<MatrixPolicy>*>(state.temporary_variables_.get());
+        static_cast<RosenbrockTemporaryVariables<DenseMatrixPolicy>*>(state.temporary_variables_.get());
     auto& Ynew = derived_class_temporary_variables->Ynew_;
     auto& initial_forcing = derived_class_temporary_variables->initial_forcing_;
     auto& K = derived_class_temporary_variables->K_;
@@ -71,14 +72,19 @@ namespace micm
       //  Repeat step calculation until current step accepted
       while (!accepted)
       {
-        // Compute alpha accounting for the last alpha value
-        // This is necessary to avoid the need to re-factor the jacobian
-        double alpha = 1.0 / (H * parameters.gamma_[0]) - last_alpha;
+        // Compute alpha for AlphaMinusJacobian function 
+        double alpha = 1.0 / (H * parameters.gamma_[0]);
+        if constexpr (!LinearSolverInPlaceConcept<LinearSolverPolicy, DenseMatrixPolicy, SparseMatrixPolicy>)
+        { 
+          // Compute alpha accounting for the last alpha value
+          // This is necessary to avoid the need to re-factor the jacobian for non-inline LU algorithms
+          alpha -= last_alpha;
+          last_alpha = alpha;
+        }
 
         // Form and factor the rosenbrock ode jacobian
-        LinearFactor(alpha, Y, stats, state);
-        last_alpha = alpha;
-
+        LinearFactor(alpha, stats, state);
+
         // Compute the stages
         for (uint64_t stage = 0; stage < parameters.stages_; ++stage)
         {
@@ -109,7 +115,14 @@ namespace micm
           {
             K[stage].Axpy(parameters.c_[stage_combinations + j] / H, K[j]);
           }
-          linear_solver_.Solve(K[stage], state.lower_matrix_, state.upper_matrix_);
+          if constexpr (LinearSolverInPlaceConcept<LinearSolverPolicy, DenseMatrixPolicy, SparseMatrixPolicy>)
+          {
+            linear_solver_.Solve(K[stage], state.jacobian_);
+          }
+          else
+          {
+            linear_solver_.Solve(K[stage], state.lower_matrix_, state.upper_matrix_);
+          }
           stats.solves_ += 1;
         }
 
@@ -178,6 +191,13 @@ namespace micm
           {
             stats.rejected_ += 1;
           }
+          // Re-generate the Jacobian matrix for the inline LU algorithm
+          if constexpr (LinearSolverInPlaceConcept<LinearSolverPolicy, DenseMatrixPolicy, SparseMatrixPolicy>)
+          {
+            state.jacobian_.Fill(0);
+            rates_.SubtractJacobianTerms(state.rate_constants_, Y, state.jacobian_);
+            stats.jacobian_updates_ += 1; 
+          }
         }
       }
     }
@@ -232,15 +252,23 @@ namespace micm
   template<class RatesPolicy, class LinearSolverPolicy, class Derived>
   inline void AbstractRosenbrockSolver<RatesPolicy, LinearSolverPolicy, Derived>::LinearFactor(
       const double alpha,
-      const auto& number_densities,
       SolverStats& stats,
       auto& state) const
   {
     MICM_PROFILE_FUNCTION();
-
+    using DenseMatrixPolicy = decltype(state.variables_);
+    using SparseMatrixPolicy = decltype(state.jacobian_);
+
     static_cast<const Derived*>(this)->AlphaMinusJacobian(state.jacobian_, alpha);
 
-    linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_);
+    if constexpr (LinearSolverInPlaceConcept<LinearSolverPolicy, DenseMatrixPolicy, SparseMatrixPolicy>)
+    {
+      linear_solver_.Factor(state.jacobian_);
+    }
+    else
+    {
+      linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_);
+    }
     stats.decompositions_ += 1;
   }