diff --git a/include/micm/solver/solver.hpp b/include/micm/solver/solver.hpp index 5b3a95d8a..7a8f9fb7e 100644 --- a/include/micm/solver/solver.hpp +++ b/include/micm/solver/solver.hpp @@ -13,13 +13,17 @@ namespace micm { - class SolverImplBase { - public: + using SolverParameters = std::variant; + + class SolverImplBase + { + public: virtual ~SolverImplBase() = default; }; template - struct SolverImpl : public SolverImplBase { + struct SolverImpl : public SolverImplBase + { ProcessSetPolicy process_set_; LinearSolverPolicy linear_solver_; }; @@ -27,9 +31,6 @@ namespace micm class Solver { private: - using SolverParameters = - std::variant; - SolverParameters parameters_; std::size_t number_of_grid_cells_; std::size_t number_of_species_; diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp index 5e0694257..2d306f0f3 100644 --- a/include/micm/solver/solver_builder.hpp +++ b/include/micm/solver/solver_builder.hpp @@ -26,70 +26,68 @@ namespace micm micm::System system_; std::size_t number_of_grid_cells_; std::vector reactions_; - std::variant options_; + SolverParameters options_; bool ignore_unused_species_ = true; bool reorder_state_ = true; public: /// @brief Set the chemical system - /// @param system - /// @return + /// @param system + /// @return SolverBuilder& SetSystem(micm::System system); /// @brief Set the reactions - /// @param reactions - /// @return + /// @param reactions + /// @return SolverBuilder& SetReactions(std::vector reactions); /// @brief Set the number of grid cells - /// @param number_of_grid_cells - /// @return + /// @param number_of_grid_cells + /// @return SolverBuilder& SetNumberOfGridCells(int number_of_grid_cells); /// @brief Choose a rosenbrock solver - /// @param options - /// @return - SolverBuilder& SolverParameters(const RosenbrockSolverParameters& options); + /// @param options + /// @return + SolverBuilder& SetSolverParameters(const RosenbrockSolverParameters& options); /// @brief Choose a backward euler solver - /// @param options - /// @return - SolverBuilder& SolverParameters(const BackwardEulerSolverParameters& options); + /// @param options + /// @return + SolverBuilder& SetSolverParameters(const BackwardEulerSolverParameters& options); - /// @brief - /// @return + /// @brief + /// @return Solver Build(); protected: - - /// @brief - /// @return + /// @brief + /// @return virtual Solver BuildBackwardEulerSolver() = 0; virtual ~SolverBuilder() = default; - /// @brief - /// @tparam ProcessSetPolicy - /// @return + /// @brief + /// @tparam ProcessSetPolicy + /// @return template void UnusedSpeciesCheck(); /// @brief Get a species map properly ordered - /// @return + /// @return template std::map GetSpeciesMap() const; /// @brief Set the absolute tolerances per species - /// @param parameters - /// @param species_map - /// @return + /// @param parameters + /// @param species_map + /// @return void SetAbsoluteTolerances(std::vector& tolerances, const std::map& species_map) const; /// @brief Return the labels of the custom parameters - /// @return + /// @return std::vector GetCustomParameterLabels() const; std::vector GetJacobianDiagonalElements(auto jacobian) const; - }; template diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index c9797885f..e030f9e55 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -65,7 +65,7 @@ namespace micm return *this; } - inline SolverBuilder& SolverBuilder::SolverParameters(const RosenbrockSolverParameters& options) + inline SolverBuilder& SolverBuilder::SetSolverParameters(const RosenbrockSolverParameters& options) { if (!std::holds_alternative(options_)) { @@ -75,7 +75,7 @@ namespace micm return *this; } - inline SolverBuilder& SolverBuilder::SolverParameters(const BackwardEulerSolverParameters& options) + inline SolverBuilder& SolverBuilder::SetSolverParameters(const BackwardEulerSolverParameters& options) { if (!std::holds_alternative(options_)) { @@ -137,7 +137,6 @@ namespace micm for (auto& name : system_.UniqueNames()) species_map[name] = index++; - if (reorder_state_) { // get unsorted Jacobian non-zero elements @@ -189,7 +188,8 @@ namespace micm } } - inline std::vector SolverBuilder::GetCustomParameterLabels() const { + inline std::vector SolverBuilder::GetCustomParameterLabels() const + { std::vector param_labels{}; for (const auto& reaction : reactions_) if (reaction.rate_constant_) @@ -198,7 +198,8 @@ namespace micm return param_labels; } - inline std::vector SolverBuilder::GetJacobianDiagonalElements(auto jacobian) const { + inline std::vector SolverBuilder::GetJacobianDiagonalElements(auto jacobian) const + { std::vector jacobian_diagonal_elements; jacobian_diagonal_elements.reserve(jacobian.NumRows()); @@ -221,19 +222,23 @@ namespace micm auto labels = GetCustomParameterLabels(); std::size_t number_of_species = system_.StateSize(); - UnusedSpeciesCheck(); SetAbsoluteTolerances(parameters.absolute_tolerance_, species_map); ProcessSetPolicy process_set(reactions_, species_map); - auto jacobian = BuildJacobian(process_set.NonZeroJacobianElements(), number_of_grid_cells_, number_of_species); + auto jacobian = + BuildJacobian(process_set.NonZeroJacobianElements(), number_of_grid_cells_, number_of_species); auto diagonal_elements = GetJacobianDiagonalElements(jacobian); + process_set.SetJacobianFlatIds(jacobian); - micm::LinearSolver linear_solver(jacobian, 1e-30); + micm::LinearSolver linear_solver(jacobian, 1e-30); return Solver( - new SolverImpl(), - parameters, number_of_grid_cells_, number_of_species, reactions_.size()); + new SolverImpl(), + parameters, + number_of_grid_cells_, + number_of_species, + reactions_.size()); } } // namespace micm \ No newline at end of file diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp index affc90ed5..d47baf727 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -9,7 +9,8 @@ #include #include -namespace { +namespace +{ auto a = micm::Species("A"); auto b = micm::Species("B"); auto c = micm::Species("C"); @@ -30,7 +31,7 @@ namespace { .SetPhase(gas_phase); micm::System the_system = micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }); std::vector reactions = { r1, r2 }; -} +} // namespace TEST(SolverBuilder, CanBuildBackwardEuler) { @@ -38,16 +39,16 @@ TEST(SolverBuilder, CanBuildBackwardEuler) .SetSystem(the_system) .SetReactions(reactions) .SetNumberOfGridCells(1) - .SolverParameters(micm::BackwardEulerSolverParameters{}) + .SetSolverParameters(micm::BackwardEulerSolverParameters{}) .Build(); constexpr std::size_t L = 4; auto backward_euler_vector = micm::CpuSolverBuilder, micm::SparseMatrix>() - .SetSystem(the_system) - .SetReactions(reactions) - .SetNumberOfGridCells(1) - .SolverParameters(micm::BackwardEulerSolverParameters{}) - .Build(); + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(1) + .SetSolverParameters(micm::BackwardEulerSolverParameters{}) + .Build(); } TEST(SolverBuilder, CanBuildRosenbrock)