hiding some template details behind a solver impl

include/micm/solver/solver.hpp

@@ -8,10 +8,22 @@
 #include <micm/solver/backward_euler_solver_parameters.hpp>
 #include <micm/solver/rosenbrock_solver_parameters.hpp>
 
+#include <memory>
 #include <variant>
 
 namespace micm
 {
+  class SolverImplBase {
+  public:
+    virtual ~SolverImplBase() = default;
+  };
+
+  template<class LinearSolverPolicy, class ProcessSetPolicy>
+  struct SolverImpl : public SolverImplBase {
+    ProcessSetPolicy process_set_;
+    LinearSolverPolicy linear_solver_;
+  };
+
   class Solver
   {
    private:
@@ -22,6 +34,7 @@ namespace micm
     std::size_t number_of_grid_cells_;
     std::size_t number_of_species_;
     std::size_t number_of_reactions_;
+    std::unique_ptr<SolverImplBase> solver_impl_;
 
    public:
     Solver()
@@ -33,11 +46,13 @@ namespace micm
     }
 
     Solver(
+        SolverImplBase* solver_impl,
         SolverParameters parameters,
         std::size_t number_of_grid_cells,
         std::size_t number_of_species,
         std::size_t number_of_reactions)
-        : parameters_(parameters),
+        : solver_impl_(solver_impl),
+          parameters_(parameters),
           number_of_grid_cells_(number_of_grid_cells),
           number_of_species_(number_of_species),
           number_of_reactions_(number_of_reactions)

include/micm/solver/solver_builder.hpp

@@ -85,6 +85,12 @@ namespace micm
     /// @return 
     void SetAbsoluteTolerances(std::vector<double>& tolerances, const std::map<std::string, std::size_t>& species_map) const;
 
+    /// @brief Return the labels of the custom parameters
+    /// @return 
+    std::vector<std::string> GetCustomParameterLabels() const;
+
+    std::vector<std::size_t> GetJacobianDiagonalElements(auto jacobian) const;
+
   };
 }  // namespace micm
 

include/micm/solver/solver_builder.inl

@@ -104,15 +104,22 @@ namespace micm
   Solver SolverBuilder::BuildBackwardEulerSolver()
   {
     auto parameters = std::get<BackwardEulerSolverParameters>(options_);
-    UnusedSpeciesCheck<micm::ProcessSet>();
     auto species_map = GetSpeciesMap<MatrixPolicy, micm::ProcessSet>();
-    SetAbsoluteTolerances(parameters.absolute_tolerance_, species_map);
+    auto labels = GetCustomParameterLabels();
     std::size_t number_of_species = system_.StateSize();
 
+    UnusedSpeciesCheck<micm::ProcessSet>();
+    SetAbsoluteTolerances(parameters.absolute_tolerance_, species_map);
+
     micm::ProcessSet process_set(reactions_, species_map);
     auto jacobian = BuildJacobian<SparseMatrixPolicy>(process_set.NonZeroJacobianElements(), number_of_grid_cells_, number_of_species);
+    auto diagonal_elements = GetJacobianDiagonalElements(jacobian);
+    process_set.SetJacobianFlatIds(jacobian);
+    micm::LinearSolver<double, SparseMatrixPolicy> linear_solver(jacobian, 1e-30);
 
-    return Solver(parameters, number_of_grid_cells_, number_of_species, reactions_.size());
+    return Solver(
+      new SolverImpl<decltype(linear_solver), decltype(process_set)>(), 
+      parameters, number_of_grid_cells_, number_of_species, reactions_.size());
   }
 
   template<class ProcessSetPolicy>
@@ -203,4 +210,26 @@ namespace micm
       }
     }
   }
+
+  std::vector<std::string> SolverBuilder::GetCustomParameterLabels() const {
+    std::vector<std::string> param_labels{};
+    for (const auto& reaction : reactions_)
+      if (reaction.rate_constant_)
+        for (auto& label : reaction.rate_constant_->CustomParameters())
+          param_labels.push_back(label);
+    return param_labels;
+  }
+
+  std::vector<std::size_t> SolverBuilder::GetJacobianDiagonalElements(auto jacobian) const {
+    std::vector<std::size_t> jacobian_diagonal_elements;
+
+    jacobian_diagonal_elements.reserve(jacobian.NumRows());
+
+    for (std::size_t i = 0; i < jacobian.NumRows(); ++i)
+    {
+      jacobian_diagonal_elements.push_back(jacobian.VectorIndex(0, i, i));
+    }
+
+    return jacobian_diagonal_elements;
+  }
 }  // namespace micm