removed compile errors

NCAR · Nov 7, 2024 · 1b67f79 · 1b67f79
1 parent 4c84691
commit 1b67f79
Show file tree

Hide file tree

Showing 6 changed files with 50 additions and 51 deletions.
diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp
@@ -49,7 +49,8 @@ namespace micm
         const BackwardEulerSolverParameters& parameters,
         LinearSolverPolicy&& linear_solver,
         RatesPolicy&& rates,
-        auto& jacobian)
+        auto& jacobian,
+        const size_t number_of_species)
         : parameters_(parameters),
           linear_solver_(std::move(linear_solver)),
           rates_(std::move(rates)),
@@ -78,14 +79,14 @@ namespace micm
     /// @return true if the residual is small enough to stop the iteration
     template<class DenseMatrixPolicy>
     static bool IsConverged(
-        const BackwardEulerSolverParameters& parameters,
-        const DenseMatrixPolicy& residual,
-        const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense<DenseMatrixPolicy>);
+      const BackwardEulerSolverParameters& parameters,
+      const DenseMatrixPolicy& residual,
+      const DenseMatrixPolicy& Yn1, std::vector<double>& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense<DenseMatrixPolicy>);
     template<class DenseMatrixPolicy>
     static bool IsConverged(
-        const BackwardEulerSolverParameters& parameters,
-        const DenseMatrixPolicy& residual,
-        const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense<DenseMatrixPolicy>);
+      const BackwardEulerSolverParameters& parameters,
+      const DenseMatrixPolicy& residual,
+      const DenseMatrixPolicy& Yn1, std::vector<double>& absolute_tolerance, double relative_tolerance) requires(VectorizableDense<DenseMatrixPolicy>);
   };
 
 }  // namespace micm

diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl
@@ -133,7 +133,7 @@ namespace micm
           continue;
 
         // check for convergence
-        converged = IsConverged(parameters_, forcing, Yn1, state);
+        converged = IsConverged(parameters_, forcing, Yn1, state.absolute_tolerance_.AsVector(), state.relative_tolerance_);
       } while (!converged && iterations < max_iter);
 
       if (!converged)
@@ -180,23 +180,23 @@ namespace micm
   inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(
       const BackwardEulerSolverParameters& parameters,
       const DenseMatrixPolicy& residual,
-      const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& Yn1, std::vector<double>& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense<DenseMatrixPolicy>)
   {
     double small = parameters.small_;
-    double rel_tol = stateParams.relative_tolerance_;
-    auto& abs_tol = stateParams.absolute_tolerance_;
+    double rel_tol = relative_tolerance;
+    auto& abs_tol = absolute_tolerance;
     auto residual_iter = residual.AsVector().begin();
-    auto state_iter = state.AsVector().begin();
+    auto Yn1_iter = Yn1.AsVector().begin();
     const std::size_t n_elem = residual.NumRows() * residual.NumColumns();
     const std::size_t n_vars = abs_tol.size();
     for (std::size_t i = 0; i < n_elem; ++i)
     {
       if (std::abs(*residual_iter) > small && std::abs(*residual_iter) > abs_tol[i % n_vars] &&
-          std::abs(*residual_iter) > rel_tol * std::abs(*state_iter))
+          std::abs(*residual_iter) > rel_tol * std::abs(*Yn1_iter))
       {
         return false;
       }
-      ++residual_iter, ++state_iter;
+      ++residual_iter, ++Yn1_iter;
     }
     return true;
   }
@@ -206,13 +206,13 @@ namespace micm
   inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(
       const BackwardEulerSolverParameters& parameters,
       const DenseMatrixPolicy& residual,
-      const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& Yn1, std::vector<double>& absolute_tolerance, double relative_tolerance) requires(VectorizableDense<DenseMatrixPolicy>)
   {
     double small = parameters.small_;
-    double rel_tol = stateParams.relative_tolerance_;
-    auto& abs_tol = stateParams.absolute_tolerance_;
+    double rel_tol = relative_tolerance;
+    auto& abs_tol = absolute_tolerance;
     auto residual_iter = residual.AsVector().begin();
-    auto state_iter = state.AsVector().begin();
+    auto Yn1_iter = Yn1.AsVector().begin();
     const std::size_t n_elem = residual.NumRows() * residual.NumColumns();
     const std::size_t L = residual.GroupVectorSize();
     const std::size_t n_vars = abs_tol.size();
@@ -221,11 +221,11 @@ namespace micm
     for (std::size_t i = 0; i < whole_blocks; ++i)
     {
       if (std::abs(*residual_iter) > small && std::abs(*residual_iter) > abs_tol[(i / L) % n_vars] &&
-          std::abs(*residual_iter) > rel_tol * std::abs(*state_iter))
+          std::abs(*residual_iter) > rel_tol * std::abs(*Yn1_iter))
       {
         return false;
       }
-      ++residual_iter, ++state_iter;
+      ++residual_iter, ++Yn1_iter;
     }
 
     // evaluate the remaining rows
@@ -238,7 +238,7 @@ namespace micm
         for (std::size_t i = offset; i < offset + remaining_rows; ++i)
         {
           if (std::abs(residual_iter[i]) > small && std::abs(residual_iter[i]) > abs_tol[y] &&
-              std::abs(residual_iter[i]) > rel_tol * std::abs(state_iter[i]))
+              std::abs(residual_iter[i]) > rel_tol * std::abs(Yn1_iter[i]))
           {
             return false;
           }

diff --git a/include/micm/solver/backward_euler_solver_parameters.hpp b/include/micm/solver/backward_euler_solver_parameters.hpp
@@ -19,8 +19,6 @@ namespace micm
     template<class RatesPolicy, class LinearSolverPolicy>
     using SolverType = BackwardEuler<RatesPolicy, LinearSolverPolicy>;
 
-    std::vector<double> absolute_tolerance_;
-    double relative_tolerance_{ 1.0e-8 };
     double small_{ 1.0e-40 };
     size_t max_number_of_steps_{ 11 };
     // The time step reductions are used to determine the time step after a failed solve

diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp
@@ -1332,7 +1332,7 @@ void test_analytical_robertson(
 
   auto state = solver.GetState();
   state.relative_tolerance_ = 1e-10;
-  state.absolute_tolerance_ = std::vector<double>(3, state.relative_tolerance_ * 1e-2);
+  state.absolute_tolerance_.AsVector() = std::vector<double>(3, state.relative_tolerance_ * 1e-2);
 
   double k1 = 0.04;
   double k2 = 3e7;
@@ -1548,7 +1548,7 @@ void test_analytical_oregonator(
   auto state = solver.GetState();
 
   state.relative_tolerance_ = 1e-6;
-  state.absolute_tolerance_ = std::vector<double>(5, state.relative_tolerance_ * 1e-2);
+  state.absolute_tolerance_.AsVector() = std::vector<double>(5, state.relative_tolerance_ * 1e-2);
 
   state.SetCustomRateParameter("r1", 1.34 * 0.06);
   state.SetCustomRateParameter("r2", 1.6e9);
@@ -1726,7 +1726,7 @@ void test_analytical_hires(
 
   auto state = solver.GetState();
   state.relative_tolerance_ = 1e-6;
-  state.absolute_tolerance_ = std::vector<double>(8, state.relative_tolerance_ * 1e-2);
+  state.absolute_tolerance_.AsVector() = std::vector<double>(8, state.relative_tolerance_ * 1e-2);
 
   state.SetCustomRateParameter("r1", 1.71);
   state.SetCustomRateParameter("r2", 8.75);
@@ -1886,10 +1886,10 @@ void test_analytical_e5(
   auto state = solver.GetState();
 
   state.relative_tolerance_ = 1e-13;
-  state.absolute_tolerance_ = std::vector<double>(6, 1e-17);
-  state.absolute_tolerance_[0] = 1e-7;
-  state.absolute_tolerance_[4] = 1e-7;
-  state.absolute_tolerance_[5] = 1e-7;
+  state.absolute_tolerance_.AsVector() = std::vector<double>(6, 1e-17);
+  state.absolute_tolerance_[0][0] = 1e-7;
+  state.absolute_tolerance_[0][4] = 1e-7;
+  state.absolute_tolerance_[0][5] = 1e-7;
 
   state.SetCustomRateParameter("r1", 7.89e-10);
   state.SetCustomRateParameter("r2", 1.13e9);

diff --git a/test/integration/test_chapman_integration.cpp b/test/integration/test_chapman_integration.cpp
@@ -43,17 +43,18 @@ TEST(ChapmanIntegration, CanBuildChapmanSystemUsingConfig)
   micm::State state = solver.GetState();
   state.SetRelativeTolerances(solver_params.relative_tolerance_); 
   EXPECT_EQ(state.relative_tolerance_, 1.0e-4);
+  auto& abs_tol = state.absolute_tolerance_.AsVector();
 
   for (size_t n_grid_cell = 0; n_grid_cell < state.number_of_grid_cells_; ++n_grid_cell)
   {
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["Ar"]], 1.0e-12);
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["CO2"]], 1.0e-12);
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["H2O"]], 1.0e-12);
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["N2"]], 1.0e-12);
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O1D"]], 1.0e-12);
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O"]], 1.0e-12);
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O2"]], 1.0e-12);
-    EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O3"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["Ar"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["CO2"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["H2O"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["N2"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["O1D"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["O"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["O2"]], 1.0e-12);
+    EXPECT_EQ(abs_tol[state.variable_map_["O3"]], 1.0e-12);
   }
 
   // User gives an input of concentrations

diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp
@@ -68,26 +68,25 @@ void CheckIsConverged()
 
   micm::BackwardEulerSolverParameters parameters;
   DenseMatrixPolicy residual{ 4, 3, 0.0 };
-  DenseMatrixPolicy state{ 4, 3, 0.0 };
-  micm::StateParameters state_params;
+  DenseMatrixPolicy Yn1{ 4, 3, 0.0 };
 
   parameters.small_ = 1e-6;
-  state_params.relative_tolerance_ = 1e-3;
-  state_params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 };
+  double relative_tolerance = 1e-3;
+  std::vector<double> absolute_tolerance = { 1e-6, 1e-6, 1e-6 };
 
-  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params));
+  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance));
   residual[0][1] = 1e-5;
-  ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params));
+  ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance));
   parameters.small_ = 1e-4;
-  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params));
+  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance));
   residual[3][2] = 1e-3;
-  ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params));
-  state[3][2] = 10.0;
-  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params));
+  ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance));
+  Yn1[3][2] = 10.0;
+  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance));
   residual[3][2] = 1e-1;
-  ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params));
-  state_params.absolute_tolerance_[2] = 1.0;
-  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params));  
+  ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance));
+  absolute_tolerance[2] = 1.0;
+  ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance));  
 }
 
 TEST(BackwardEuler, IsConverged)