From 3455e945249ccd86a40fccfd6a058a7b253aa5fa Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Wed, 5 Jun 2024 12:19:20 -0700 Subject: [PATCH] apply solver builder to jit solver --- include/micm/solver/backward_euler.hpp | 17 +-- include/micm/solver/jit_rosenbrock.hpp | 75 ++++++------- include/micm/solver/jit_solver_builder.hpp | 11 +- include/micm/solver/jit_solver_parameters.hpp | 25 +++++ include/micm/solver/linear_solver.hpp | 5 + include/micm/solver/lu_decomposition.hpp | 21 ++++ include/micm/solver/lu_decomposition.inl | 5 + include/micm/solver/rosenbrock.hpp | 13 ++- include/micm/solver/solver.hpp | 4 +- include/micm/solver/solver_builder.hpp | 4 - include/micm/solver/solver_builder.inl | 20 +--- include/micm/solver/state.hpp | 1 - include/micm/util/sparse_matrix.hpp | 13 +++ .../integration/analytical_jit_rosenbrock.cpp | 57 +++++----- test/integration/jit_terminator.cpp | 44 ++++---- .../regression/RosenbrockChapman/jit_util.hpp | 86 ++++++++------- .../regression_test_jit_dforce_dy.cpp | 6 +- .../regression_test_jit_p_force.cpp | 15 +-- .../regression_test_jit_solve.cpp | 30 +----- test/tutorial/test_jit_tutorial.cpp | 18 ++-- test/unit/openmp/run_solver.hpp | 1 - test/unit/openmp/test_openmp.cpp | 6 +- test/unit/openmp/test_openmp_jit_solver.cpp | 17 ++- test/unit/solver/test_jit_rosenbrock.cpp | 88 +++++---------- test/unit/solver/test_solver_builder.cpp | 15 +-- test/unit/util/test_sparse_matrix_policy.hpp | 101 +++++++++++------- .../test_sparse_matrix_standard_ordering.cpp | 24 ++--- .../test_sparse_matrix_vector_ordering.cpp | 30 +++--- 28 files changed, 399 insertions(+), 353 deletions(-) create mode 100644 include/micm/solver/jit_solver_parameters.hpp diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index 7f894b543..94d87105d 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/include/micm/solver/backward_euler.hpp @@ -46,18 +46,23 @@ namespace micm /// @brief Default constructor BackwardEuler( BackwardEulerSolverParameters parameters, - LinearSolverPolicy linear_solver, - ProcessSetPolicy process_set, - std::vector jacobian_diagonal_elements, + LinearSolverPolicy&& linear_solver, + ProcessSetPolicy&& process_set, + auto& jacobian, std::vector& processes) : parameters_(parameters), - linear_solver_(linear_solver), - process_set_(process_set), - jacobian_diagonal_elements_(jacobian_diagonal_elements), + linear_solver_(std::move(linear_solver)), + process_set_(std::move(process_set)), + jacobian_diagonal_elements_(jacobian.DiagonalIndices(0)), processes_(processes) { } + BackwardEuler(const BackwardEuler&) = delete; + BackwardEuler& operator=(const BackwardEuler&) = delete; + BackwardEuler(BackwardEuler&&) = default; + BackwardEuler& operator=(BackwardEuler&&) = default; + virtual ~BackwardEuler() = default; /// @brief Advances the given step over the specified time step diff --git a/include/micm/solver/jit_rosenbrock.hpp b/include/micm/solver/jit_rosenbrock.hpp index 9f0a1ea0d..e2275f0d8 100644 --- a/include/micm/solver/jit_rosenbrock.hpp +++ b/include/micm/solver/jit_rosenbrock.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -30,22 +31,24 @@ namespace micm { - template< - template class MatrixPolicy = VectorMatrix, - class SparseMatrixPolicy = DefaultVectorSparseMatrix, - class LinearSolverPolicy = JitLinearSolver, - class ProcessSetPolicy = JitProcessSet> - class JitRosenbrockSolver : public RosenbrockSolver + /// @brief A Rosenbrock solver with JIT-compiled optimizations + + template + class JitRosenbrockSolver : public RosenbrockSolver { llvm::orc::ResourceTrackerSP function_resource_tracker_; using FuncPtr = void (*)(double*, const double); FuncPtr alpha_minus_jacobian_ = nullptr; public: + + /// @brief Solver parameters typename + using ParametersType = RosenbrockSolverParameters; + JitRosenbrockSolver(const JitRosenbrockSolver&) = delete; JitRosenbrockSolver& operator=(const JitRosenbrockSolver&) = delete; JitRosenbrockSolver(JitRosenbrockSolver&& other) - : RosenbrockSolver(std::move(other)), + : RosenbrockSolver(std::move(other)), function_resource_tracker_(std::move(other.function_resource_tracker_)), alpha_minus_jacobian_(std::move(other.alpha_minus_jacobian_)) { @@ -54,7 +57,7 @@ namespace micm JitRosenbrockSolver& operator=(JitRosenbrockSolver&& other) { - RosenbrockSolver::operator=(std::move(other)); + RosenbrockSolver::operator=(std::move(other)); function_resource_tracker_ = std::move(other.function_resource_tracker_); alpha_minus_jacobian_ = std::move(other.alpha_minus_jacobian_); other.alpha_minus_jacobian_ = NULL; @@ -62,34 +65,25 @@ namespace micm } /// @brief Builds a Rosenbrock solver for the given system, processes, and solver parameters - /// @param system The chemical system to create the solver for - /// @param processes The collection of chemical processes that will be applied during solving + /// @param parameters Solver parameters + /// @param linear_solver Linear solver + /// @param process_set Process set + /// @param jacobian Jacobian matrix + /// @param processes Vector of processes JitRosenbrockSolver( - const System& system, - const std::vector& processes, - const RosenbrockSolverParameters& parameters) - : RosenbrockSolver( - system, - processes, + RosenbrockSolverParameters parameters, + LinearSolverPolicy linear_solver, + ProcessSetPolicy process_set, + auto& jacobian, + std::vector& processes) + : RosenbrockSolver( parameters, - [&](const SparseMatrixPolicy& matrix, double initial_value) -> LinearSolverPolicy { - return LinearSolverPolicy{ matrix, initial_value }; - }, - [&](const std::vector& processes, - const std::map& variable_map) -> ProcessSetPolicy { - return ProcessSetPolicy{ processes, variable_map }; - }) + std::move(linear_solver), + std::move(process_set), + jacobian, + processes) { - MatrixPolicy temp{}; - if (temp.GroupVectorSize() != parameters.number_of_grid_cells_) - { - std::string msg = - "JIT functions require the number of grid cells solved together to match the vector dimension template " - "parameter, currently: " + - std::to_string(temp.GroupVectorSize()); - throw std::system_error(make_error_code(MicmJitErrc::InvalidMatrix), msg); - } - this->GenerateAlphaMinusJacobian(); + this->GenerateAlphaMinusJacobian(jacobian); } ~JitRosenbrockSolver() @@ -104,8 +98,17 @@ namespace micm /// @brief compute [alpha * I - dforce_dy] /// @param jacobian Jacobian matrix (dforce_dy) /// @param alpha + template void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const { + if (jacobian.GroupVectorSize() != jacobian.NumberOfBlocks()) + { + std::string msg = + "JIT functions require the number of grid cells solved together to match the vector dimension template " + "parameter, currently: " + + std::to_string(jacobian.GroupVectorSize()); + throw std::system_error(make_error_code(MicmJitErrc::InvalidMatrix), msg); + } double a = alpha; if (alpha_minus_jacobian_) { @@ -119,9 +122,9 @@ namespace micm } private: - void GenerateAlphaMinusJacobian() + void GenerateAlphaMinusJacobian(auto& jacobian) { - auto jacobian = this->GetState().jacobian_; + auto diagonal_elements = jacobian.DiagonalIndices(0); // save sizes needed throughout the function std::size_t n_cells = jacobian.GroupVectorSize(); std::size_t number_of_nonzero_jacobian_elements = jacobian.AsVector().size(); @@ -143,7 +146,7 @@ namespace micm // iterative over the blocks of the jacobian and add the alpha value // jacobian_vector[i_elem + i_cell] += alpha; - for (const auto& i_elem : this->state_parameters_.jacobian_diagonal_elements_) + for (const auto& i_elem : diagonal_elements) { llvm::Value* ptr_index[1]; diff --git a/include/micm/solver/jit_solver_builder.hpp b/include/micm/solver/jit_solver_builder.hpp index cdff5ed8a..bc1ae16ff 100644 --- a/include/micm/solver/jit_solver_builder.hpp +++ b/include/micm/solver/jit_solver_builder.hpp @@ -11,10 +11,15 @@ namespace micm { /// @brief Builder of CPU-based solvers optimized for a particular chemical system using JIT compilation /// @tparam SolverParametersPolicy Policy for the ODE solver - /// @tparam T Primary data type /// @tparam L Vector size /// /// JIT-compiled solvers only work with vector-ordered matrices - template - using JitSolverBuilder = SolverBuilder, SparseMatrix>, JitProcessSet, JitLinearSolver>, JitLuDecomposition>>; + template + using JitSolverBuilder = SolverBuilder, + SparseMatrix>, + JitProcessSet, + JitLinearSolver>, + JitLuDecomposition>>; } // namespace micm \ No newline at end of file diff --git a/include/micm/solver/jit_solver_parameters.hpp b/include/micm/solver/jit_solver_parameters.hpp new file mode 100644 index 000000000..06ae8db99 --- /dev/null +++ b/include/micm/solver/jit_solver_parameters.hpp @@ -0,0 +1,25 @@ +// Copyright (C) 2023-2024 National Center for Atmospheric Research +// SPDX-License-Identifier: Apache-2.0 +#pragma once + +#include "rosenbrock_solver_parameters.hpp" +#include +#include + +namespace micm +{ + /// @brief Parameters for the JIT Rosenbrock solver + struct JitRosenbrockSolverParameters : public RosenbrockSolverParameters + { + template + using SolverType = JitRosenbrockSolver; + + /// @brief Constructor from base class + /// @param base + JitRosenbrockSolverParameters(const RosenbrockSolverParameters& base) + : RosenbrockSolverParameters(base) + { + } + + }; +} \ No newline at end of file diff --git a/include/micm/solver/linear_solver.hpp b/include/micm/solver/linear_solver.hpp index 6bf5e9aaa..139ddeeb7 100644 --- a/include/micm/solver/linear_solver.hpp +++ b/include/micm/solver/linear_solver.hpp @@ -57,6 +57,11 @@ namespace micm /// @brief default constructor LinearSolver(){}; + LinearSolver(const LinearSolver&) = delete; + LinearSolver& operator=(const LinearSolver&) = delete; + LinearSolver(LinearSolver&&) = default; + LinearSolver& operator=(LinearSolver&&) = default; + /// @brief Constructs a linear solver for the sparsity structure of the given matrix /// @param matrix Sparse matrix /// @param initial_value Initial value for matrix elements diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp index 538818a79..c68cddd3b 100644 --- a/include/micm/solver/lu_decomposition.hpp +++ b/include/micm/solver/lu_decomposition.hpp @@ -7,6 +7,17 @@ #include #include +namespace +{ + template + concept SparseMatrixConcept = requires(T t) + { + t.NumRows(); + t.NumColumns(); + t.NumberOfBlocks(); + }; +} + namespace micm { @@ -76,20 +87,28 @@ namespace micm /// @brief default constructor LuDecomposition(); + LuDecomposition(const LuDecomposition&) = delete; + LuDecomposition& operator=(const LuDecomposition&) = delete; + LuDecomposition(LuDecomposition&&) = default; + LuDecomposition& operator=(LuDecomposition&&) = default; + /// @brief Construct an LU decomposition algorithm for a given sparse matrix /// @param matrix Sparse matrix template + requires(SparseMatrixConcept) LuDecomposition(const SparseMatrixPolicy& matrix); /// @brief Create an LU decomposition algorithm for a given sparse matrix policy /// @param matrix Sparse matrix template + requires(SparseMatrixConcept) static LuDecomposition Create(const SparseMatrixPolicy& matrix); /// @brief Create sparse L and U matrices for a given A matrix /// @param A Sparse matrix that will be decomposed /// @return L and U Sparse matrices template + requires(SparseMatrixConcept) static std::pair GetLUMatrices( const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value); @@ -99,6 +118,7 @@ namespace micm /// @param L The lower triangular matrix created by decomposition /// @param U The upper triangular matrix created by decomposition template + requires(SparseMatrixConcept) void Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const; /// @brief Perform an LU decomposition on a given A matrix @@ -123,6 +143,7 @@ namespace micm /// @brief Initialize arrays for the LU decomposition /// @param A Sparse matrix to decompose template + requires(SparseMatrixConcept) void Initialize(const SparseMatrixPolicy& matrix, auto initial_value); }; diff --git a/include/micm/solver/lu_decomposition.inl b/include/micm/solver/lu_decomposition.inl index dc378db4f..75c77ece9 100644 --- a/include/micm/solver/lu_decomposition.inl +++ b/include/micm/solver/lu_decomposition.inl @@ -10,12 +10,14 @@ namespace micm } template + requires(SparseMatrixConcept) inline LuDecomposition::LuDecomposition(const SparseMatrixPolicy& matrix) { Initialize(matrix, typename SparseMatrixPolicy::value_type()); } template + requires(SparseMatrixConcept) inline LuDecomposition LuDecomposition::Create(const SparseMatrixPolicy& matrix) { LuDecomposition lu_decomp{}; @@ -24,6 +26,7 @@ namespace micm } template + requires(SparseMatrixConcept) inline void LuDecomposition::Initialize(const SparseMatrixPolicy& matrix, auto initial_value) { MICM_PROFILE_FUNCTION(); @@ -100,6 +103,7 @@ namespace micm } template + requires(SparseMatrixConcept) inline std::pair LuDecomposition::GetLUMatrices( const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value) @@ -162,6 +166,7 @@ namespace micm } template + requires(SparseMatrixConcept) inline void LuDecomposition::Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const { bool is_singular = false; diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index 3cc3d4df3..cad0d3c26 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -68,16 +68,21 @@ namespace micm RosenbrockSolverParameters parameters, LinearSolverPolicy linear_solver, ProcessSetPolicy process_set, - std::vector jacobian_diagonal_elements, + auto& jacobian, std::vector& processes) : parameters_(parameters), - linear_solver_(linear_solver), - process_set_(process_set), - jacobian_diagonal_elements_(jacobian_diagonal_elements), + linear_solver_(std::move(linear_solver)), + process_set_(std::move(process_set)), + jacobian_diagonal_elements_(jacobian.DiagonalIndices(0)), processes_(processes) { } + RosenbrockSolver(const RosenbrockSolver&) = delete; + RosenbrockSolver& operator=(const RosenbrockSolver&) = delete; + RosenbrockSolver(RosenbrockSolver&&) = default; + RosenbrockSolver& operator=(RosenbrockSolver&&) = default; + virtual ~RosenbrockSolver() = default; /// @brief Advances the given step over the specified time step diff --git a/include/micm/solver/solver.hpp b/include/micm/solver/solver.hpp index 2020b6c5a..9b3bd4e7b 100644 --- a/include/micm/solver/solver.hpp +++ b/include/micm/solver/solver.hpp @@ -21,12 +21,12 @@ namespace micm SolverPolicy solver_; Solver( - SolverPolicy solver, + SolverPolicy&& solver, StateParameters state_parameters, std::size_t number_of_grid_cells, std::size_t number_of_species, std::size_t number_of_reactions) - : solver_(solver), + : solver_(std::move(solver)), number_of_grid_cells_(number_of_grid_cells), number_of_species_(number_of_species), number_of_reactions_(number_of_reactions), diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp index f517df8e7..b42d70ab8 100644 --- a/include/micm/solver/solver_builder.hpp +++ b/include/micm/solver/solver_builder.hpp @@ -102,10 +102,6 @@ namespace micm /// @return The labels of the custom parameters std::vector GetCustomParameterLabels() const; - /// @brief Returns the diagonal elements of the jacobian matrix - /// @param jacobian - /// @return The diagonal elements of the jacobian matrix - std::vector GetJacobianDiagonalElements(auto jacobian) const; }; /// @brief Builder of CPU-based general solvers diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index 041a5fb4d..297c0aa2b 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -198,22 +198,6 @@ namespace micm return param_labels; } - template - inline std::vector SolverBuilder::GetJacobianDiagonalElements( - auto jacobian) const - { - std::vector 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; - } - template inline auto SolverBuilder::Build() { @@ -240,7 +224,6 @@ namespace micm ProcessSetPolicy process_set(this->reactions_, species_map); auto nonzero_elements = process_set.NonZeroJacobianElements(); auto jacobian = BuildJacobian(nonzero_elements, this->number_of_grid_cells_, number_of_species); - auto jacobian_diagonal_elements = this->GetJacobianDiagonalElements(jacobian); process_set.SetJacobianFlatIds(jacobian); LinearSolverPolicy linear_solver(jacobian, 1e-30); @@ -254,12 +237,11 @@ namespace micm .number_of_rate_constants_ = this->reactions_.size(), .variable_names_ = variable_names, .custom_rate_parameter_labels_ = labels, - .jacobian_diagonal_elements_ = jacobian_diagonal_elements, .nonzero_jacobian_elements_ = nonzero_elements }; return Solver>( SolverPolicy( - this->options_, linear_solver, process_set, jacobian_diagonal_elements, this->reactions_), + this->options_, std::move(linear_solver), std::move(process_set), jacobian, this->reactions_), state_parameters, this->number_of_grid_cells_, number_of_species, diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 4ff66d41d..790b0333e 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -31,7 +31,6 @@ namespace micm std::size_t number_of_rate_constants_{ 0 }; std::vector variable_names_{}; std::vector custom_rate_parameter_labels_{}; - std::vector jacobian_diagonal_elements_{}; std::set> nonzero_jacobian_elements_{}; }; diff --git a/include/micm/util/sparse_matrix.hpp b/include/micm/util/sparse_matrix.hpp index c1cca6b6a..c841682d5 100644 --- a/include/micm/util/sparse_matrix.hpp +++ b/include/micm/util/sparse_matrix.hpp @@ -207,6 +207,19 @@ namespace micm return VectorIndex(0, row, column); } + /// @brief Returns the indices of non-zero diagonal elements in a particular block + /// @param block_id Block index + /// @return Vector of indices of non-zero diagonal elements + std::vector DiagonalIndices(std::size_t block_id) const + { + std::vector indices; + indices.reserve(row_start_.size() - 1); + for (std::size_t i = 0; i < row_start_.size() - 1; ++i) + if (!IsZero(i, i)) + indices.push_back(VectorIndex(block_id, i, i)); + return indices; + } + bool IsZero(std::size_t row, std::size_t column) const { if (row >= row_start_.size() - 1 || column >= row_start_.size() - 1) diff --git a/test/integration/analytical_jit_rosenbrock.cpp b/test/integration/analytical_jit_rosenbrock.cpp index 65ae4452e..9414322eb 100644 --- a/test/integration/analytical_jit_rosenbrock.cpp +++ b/test/integration/analytical_jit_rosenbrock.cpp @@ -1,87 +1,94 @@ #include "analytical_policy.hpp" #include "analytical_surface_rxn_policy.hpp" -#include +#include +#include +#include #include -template -using DefaultVectorMatrix = micm::VectorMatrix; - -using DefaultSparseVectorMatrix = micm::SparseMatrix>; - -using DefaultJitRosenbrockSolver = micm::JitRosenbrockSolver< - DefaultVectorMatrix, - DefaultSparseVectorMatrix, - micm::JitLinearSolver<1, DefaultSparseVectorMatrix, micm::JitLuDecomposition<1>>, - micm::JitProcessSet<1>>; +constexpr std::size_t L = 1; TEST(AnalyticalExamplesJitRosenbrock, Troe) { - test_analytical_troe(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_troe(builder); } TEST(AnalyticalExamplesJitRosenbrock, TroeSuperStiffButAnalytical) { - test_analytical_stiff_troe(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_stiff_troe(builder); } TEST(AnalyticalExamplesJitRosenbrock, Photolysis) { - test_analytical_photolysis(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_photolysis(builder); } TEST(AnalyticalExamplesJitRosenbrock, PhotolysisSuperStiffButAnalytical) { - test_analytical_stiff_photolysis(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_stiff_photolysis(builder); } TEST(AnalyticalExamplesJitRosenbrock, TernaryChemicalActivation) { - test_analytical_ternary_chemical_activation(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_ternary_chemical_activation(builder); } TEST(AnalyticalExamplesJitRosenbrock, TernaryChemicalActivationSuperStiffButAnalytical) { - test_analytical_stiff_ternary_chemical_activation(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_stiff_ternary_chemical_activation(builder); } TEST(AnalyticalExamplesJitRosenbrock, Tunneling) { - test_analytical_tunneling(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_tunneling(builder); } TEST(AnalyticalExamplesJitRosenbrock, TunnelingSuperStiffButAnalytical) { - test_analytical_stiff_tunneling(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_stiff_tunneling(builder); } TEST(AnalyticalExamplesJitRosenbrock, Arrhenius) { - test_analytical_arrhenius(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_arrhenius(builder); } TEST(AnalyticalExamplesJitRosenbrock, ArrheniusSuperStiffButAnalytical) { - test_analytical_stiff_arrhenius(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_stiff_arrhenius(builder); } TEST(AnalyticalExamplesJitRosenbrock, Branched) { - test_analytical_branched(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_branched(builder); } TEST(AnalyticalExamplesJitRosenbrock, BranchedSuperStiffButAnalytical) { - test_analytical_stiff_branched(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_stiff_branched(builder); } TEST(AnalyticalExamplesJitRosenbrock, Robertson) { - test_analytical_robertson(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_robertson(builder); } TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) { - test_analytical_surface_rxn(); + auto builder = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + test_analytical_surface_rxn(builder); } diff --git a/test/integration/jit_terminator.cpp b/test/integration/jit_terminator.cpp index 180ca45ec..0864029a3 100644 --- a/test/integration/jit_terminator.cpp +++ b/test/integration/jit_terminator.cpp @@ -2,27 +2,14 @@ #include #include +#include +#include #include #include #include #include #include - -template class MatrixPolicy, class SparseMatrixPolicy> -void RunTerminatorTest() -{ - auto solver_params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(number_of_grid_cells, true); - solver_params.relative_tolerance_ = 1.0e-8; - solver_params.max_number_of_steps_ = 100000; - - TestTerminator, - micm::JitProcessSet>>(number_of_grid_cells, solver_params); -} - template using Group1VectorMatrix = micm::VectorMatrix; template @@ -39,8 +26,27 @@ using Group4SparseVectorMatrix = micm::SparseMatrix(); - RunTerminatorTest<2, Group2VectorMatrix, Group2SparseVectorMatrix>(); - RunTerminatorTest<3, Group3VectorMatrix, Group3SparseVectorMatrix>(); - RunTerminatorTest<4, Group4VectorMatrix, Group4SparseVectorMatrix>(); + auto parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); + parameters.relative_tolerance_ = 1.0e-8; + parameters.max_number_of_steps_ = 100000; + { + auto builder = micm::JitSolverBuilder(parameters) + .SetIgnoreUnusedSpecies(true); + TestTerminator(builder, 1); + } + { + auto builder = micm::JitSolverBuilder(parameters) + .SetIgnoreUnusedSpecies(true); + TestTerminator(builder, 2); + } + { + auto builder = micm::JitSolverBuilder(parameters) + .SetIgnoreUnusedSpecies(true); + TestTerminator(builder, 3); + } + { + auto builder = micm::JitSolverBuilder(parameters) + .SetIgnoreUnusedSpecies(true); + TestTerminator(builder, 4); + } } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/jit_util.hpp b/test/regression/RosenbrockChapman/jit_util.hpp index 8b1ded4ac..baad9882e 100644 --- a/test/regression/RosenbrockChapman/jit_util.hpp +++ b/test/regression/RosenbrockChapman/jit_util.hpp @@ -1,73 +1,79 @@ #include "util.hpp" #include +#include +#include +#include -template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> -micm::JitRosenbrockSolver -getTwoStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) +template +using JitBuilder = micm::JitSolverBuilder; + +template +auto getTwoStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { micm::Phase gas_phase = createGasPhase(); std::vector processes = createProcesses(gas_phase); - auto options = micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters(number_of_grid_cells); - options.ignore_unused_species_ = true; - - return micm::JitRosenbrockSolver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return JitBuilder(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()) + .SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) + .SetReactions(std::move(processes)) + .SetNumberOfGridCells(number_of_grid_cells) + .SetIgnoreUnusedSpecies(true) + .Build(); } -template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> -micm::JitRosenbrockSolver -getThreeStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) +template +auto getThreeStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { micm::Phase gas_phase = createGasPhase(); std::vector processes = createProcesses(gas_phase); - auto options = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(number_of_grid_cells); - options.ignore_unused_species_ = true; - - return micm::JitRosenbrockSolver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return JitBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) + .SetReactions(std::move(processes)) + .SetNumberOfGridCells(number_of_grid_cells) + .SetIgnoreUnusedSpecies(true) + .Build(); } -template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> -micm::JitRosenbrockSolver -getFourStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) +template +auto getFourStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { micm::Phase gas_phase = createGasPhase(); std::vector processes = createProcesses(gas_phase); - auto options = micm::RosenbrockSolverParameters::FourStageRosenbrockParameters(number_of_grid_cells); - options.ignore_unused_species_ = true; - - return micm::JitRosenbrockSolver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return JitBuilder(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()) + .SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) + .SetReactions(std::move(processes)) + .SetNumberOfGridCells(number_of_grid_cells) + .SetIgnoreUnusedSpecies(true) + .Build(); } -template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> -micm::JitRosenbrockSolver -getFourStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) +template +auto getFourStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { micm::Phase gas_phase = createGasPhase(); std::vector processes = createProcesses(gas_phase); - auto options = micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters(number_of_grid_cells); - options.ignore_unused_species_ = true; - - return micm::JitRosenbrockSolver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return JitBuilder(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()) + .SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) + .SetReactions(std::move(processes)) + .SetNumberOfGridCells(number_of_grid_cells) + .SetIgnoreUnusedSpecies(true) + .Build(); } -template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> -micm::JitRosenbrockSolver -getSixStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) +template +auto getSixStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { micm::Phase gas_phase = createGasPhase(); std::vector processes = createProcesses(gas_phase); - auto options = micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters(number_of_grid_cells); - options.ignore_unused_species_ = true; - - return micm::JitRosenbrockSolver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return JitBuilder(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()) + .SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) + .SetReactions(std::move(processes)) + .SetNumberOfGridCells(number_of_grid_cells) + .SetIgnoreUnusedSpecies(true) + .Build(); } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_jit_dforce_dy.cpp b/test/regression/RosenbrockChapman/regression_test_jit_dforce_dy.cpp index f1eb9cfcc..d9955a556 100644 --- a/test/regression/RosenbrockChapman/regression_test_jit_dforce_dy.cpp +++ b/test/regression/RosenbrockChapman/regression_test_jit_dforce_dy.cpp @@ -8,10 +8,6 @@ TEST(RegressionJitRosenbrock, VectorJacobian) { - auto solver = getThreeStageMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); + auto solver = getThreeStageMultiCellJitChapmanSolver<3>(3); testJacobian<>(solver); } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_jit_p_force.cpp b/test/regression/RosenbrockChapman/regression_test_jit_p_force.cpp index 7e88f7902..3e1889f1f 100644 --- a/test/regression/RosenbrockChapman/regression_test_jit_p_force.cpp +++ b/test/regression/RosenbrockChapman/regression_test_jit_p_force.cpp @@ -3,25 +3,18 @@ #include #include +#include #include TEST(RegressionJitRosenbrock, VectorRateConstants) { - auto solver = getThreeStageMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); + auto solver = getThreeStageMultiCellJitChapmanSolver<3>(3); testRateConstants<>(solver); } TEST(RegressionJitRosenbrock, VectorForcing) { - auto solver = getThreeStageMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); - testForcing(solver); + auto solver = getThreeStageMultiCellJitChapmanSolver<3>(3); + testForcing>(solver); } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_jit_solve.cpp b/test/regression/RosenbrockChapman/regression_test_jit_solve.cpp index 2b7376169..df6ca1e06 100644 --- a/test/regression/RosenbrockChapman/regression_test_jit_solve.cpp +++ b/test/regression/RosenbrockChapman/regression_test_jit_solve.cpp @@ -8,50 +8,30 @@ TEST(RegressionJitRosenbrock, TwoStageSolve) { - auto solver = getTwoStageMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); + auto solver = getTwoStageMultiCellJitChapmanSolver<3>(3); testSolve<>(solver, 1.0e-2); } TEST(RegressionJitRosenbrock, ThreeStageSolve) { - auto solver = getThreeStageMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); + auto solver = getThreeStageMultiCellJitChapmanSolver<3>(3); testSolve<>(solver, 1.0e-2); } TEST(RegressionJitRosenbrock, FourStageSolve) { - auto solver = getFourStageMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); + auto solver = getFourStageMultiCellJitChapmanSolver<3>(3); testSolve<>(solver, 1.0e-2); } TEST(RegressionJitRosenbrock, FourStageDASolve) { - auto solver = getFourStageDAMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); + auto solver = getFourStageDAMultiCellJitChapmanSolver<3>(3); testSolve<>(solver, 1.0e-2); } TEST(RegressionJitRosenbrock, SixStageDASolve) { - auto solver = getSixStageDAMultiCellJitChapmanSolver< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::JitLinearSolver<3, Group3SparseVectorMatrix>, - micm::JitProcessSet<3>>(3); + auto solver = getSixStageDAMultiCellJitChapmanSolver<3>(3); testSolve<>(solver, 1.0e-2); } \ No newline at end of file diff --git a/test/tutorial/test_jit_tutorial.cpp b/test/tutorial/test_jit_tutorial.cpp index 9581b0a7f..a7a93fa86 100644 --- a/test/tutorial/test_jit_tutorial.cpp +++ b/test/tutorial/test_jit_tutorial.cpp @@ -1,7 +1,9 @@ #include #include +#include #include #include +#include #include #include @@ -50,7 +52,6 @@ auto run_solver(auto& solver) auto end = std::chrono::high_resolution_clock::now(); total_solve_time += std::chrono::duration_cast(end - start); elapsed_solve_time = result.final_time_; - state.variables_ = result.result_; total_stats.function_calls_ += result.stats_.function_calls_; total_stats.jacobian_updates_ += result.stats_.jacobian_updates_; @@ -90,21 +91,20 @@ int main(const int argc, const char* argv[]) std::vector reactions{ r1, r2 }; - auto solver_parameters = RosenbrockSolverParameters::ThreeStageRosenbrockParameters(n_grid_cells); + auto solver_parameters = RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - auto solver = micm::JitSolverBuilder(solver_parameters) + auto solver = micm::CpuSolverBuilder(solver_parameters) .SetSystem(chemical_system) .SetReactions(reactions) .SetNumberOfGridCells(n_grid_cells) .Build(); auto start = std::chrono::high_resolution_clock::now(); - JitRosenbrockSolver< - GroupVectorMatrix, - GroupSparseVectorMatrix, - JitLinearSolver, - JitProcessSet> - jit_solver(chemical_system, reactions, solver_parameters); + auto jit_solver = micm::JitSolverBuilder(solver_parameters) + .SetSystem(chemical_system) + .SetReactions(reactions) + .SetNumberOfGridCells(n_grid_cells) + .Build(); auto end = std::chrono::high_resolution_clock::now(); auto jit_compile_time = std::chrono::duration_cast(end - start); diff --git a/test/unit/openmp/run_solver.hpp b/test/unit/openmp/run_solver.hpp index a206389a1..fc0c55cca 100644 --- a/test/unit/openmp/run_solver.hpp +++ b/test/unit/openmp/run_solver.hpp @@ -37,7 +37,6 @@ std::vector run_solver_on_thread_with_own_state(auto& solver, auto& stat { auto result = solver.Solve(time_step - elapsed_solve_time, state); elapsed_solve_time = result.final_time_; - state.variables_ = result.result_; } } diff --git a/test/unit/openmp/test_openmp.cpp b/test/unit/openmp/test_openmp.cpp index 6de1b5b46..01c3f04a6 100644 --- a/test/unit/openmp/test_openmp.cpp +++ b/test/unit/openmp/test_openmp.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -24,7 +25,10 @@ TEST(OpenMP, OneSolverManyStates) std::vector> results(n_threads); - RosenbrockSolver<> solver{ chemical_system, reactions, RosenbrockSolverParameters::ThreeStageRosenbrockParameters() }; + auto solver = CpuSolverBuilder(RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(chemical_system) + .SetReactions(reactions) + .Build(); #pragma omp parallel num_threads(n_threads) { diff --git a/test/unit/openmp/test_openmp_jit_solver.cpp b/test/unit/openmp/test_openmp_jit_solver.cpp index 2cc7acc3e..13f198f06 100644 --- a/test/unit/openmp/test_openmp_jit_solver.cpp +++ b/test/unit/openmp/test_openmp_jit_solver.cpp @@ -3,6 +3,8 @@ #include #include #include +#include +#include #include #include @@ -10,9 +12,8 @@ using namespace micm; -template -using Group1VectorMatrix = micm::VectorMatrix; -using Group1SparseVectorMatrix = micm::SparseMatrix>; +template +using JitBuilder = JitSolverBuilder; TEST(OpenMP, JITOneSolverManyStates) { @@ -30,12 +31,10 @@ TEST(OpenMP, JITOneSolverManyStates) std::vector> results(n_threads); - JitRosenbrockSolver< - Group1VectorMatrix, - Group1SparseVectorMatrix, - JitLinearSolver<1, Group1SparseVectorMatrix>, - micm::JitProcessSet<1>> - solver(chemical_system, reactions, RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + auto solver = JitBuilder<1>(RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(chemical_system) + .SetReactions(reactions) + .Build(); #pragma omp parallel num_threads(n_threads) { diff --git a/test/unit/solver/test_jit_rosenbrock.cpp b/test/unit/solver/test_jit_rosenbrock.cpp index b160e7887..3a34de227 100644 --- a/test/unit/solver/test_jit_rosenbrock.cpp +++ b/test/unit/solver/test_jit_rosenbrock.cpp @@ -2,6 +2,8 @@ #include #include #include +#include +#include #include #include #include @@ -10,13 +12,8 @@ #include -template class MatrixPolicy, class SparseMatrixPolicy> -micm::JitRosenbrockSolver< - MatrixPolicy, - SparseMatrixPolicy, - micm::JitLinearSolver, - micm::JitProcessSet> -getSolver() +template +SolverBuilderPolicy getSolver(SolverBuilderPolicy builder) { // ---- foo bar baz quz quuz // foo 0 1 2 - - @@ -48,24 +45,16 @@ getSolver() micm::Process r3 = micm::Process::Create().SetReactants({ quz }).SetProducts({}).SetPhase(gas_phase).SetRateConstant( micm::ArrheniusRateConstant({ .A_ = 3.5e-6 })); - return micm::JitRosenbrockSolver< - MatrixPolicy, - SparseMatrixPolicy, - micm::JitLinearSolver, - micm::JitProcessSet>( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), - std::vector{ r1, r2, r3 }, - micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(number_of_grid_cells, false)); + return builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) + .SetReactions(std::vector{ r1, r2, r3 }) + .SetReorderState(false); } -template -using SparseMatrix = micm::SparseMatrix; - -template class MatrixPolicy, class SparseMatrixPolicy> -void testAlphaMinusJacobian() +template +void testAlphaMinusJacobian(SolverBuilderPolicy builder, std::size_t number_of_grid_cells) { - auto solver = getSolver(); - // return; + builder = getSolver(builder); + auto solver = builder.SetNumberOfGridCells(number_of_grid_cells).Build(); auto jacobian = solver.GetState().jacobian_; EXPECT_EQ(jacobian.NumberOfBlocks(), number_of_grid_cells); @@ -92,7 +81,7 @@ void testAlphaMinusJacobian() jacobian[i_cell][4][2] = -53.6; jacobian[i_cell][4][4] = -1.0; } - solver.AlphaMinusJacobian(jacobian, 42.042); + solver.solver_.AlphaMinusJacobian(jacobian, 42.042); for (std::size_t i_cell = 0; i_cell < number_of_grid_cells; ++i_cell) { EXPECT_NEAR(jacobian[i_cell][0][0], 42.042 - 12.2, 1.0e-5); @@ -111,20 +100,6 @@ void testAlphaMinusJacobian() } } -template -using Group1VectorMatrix = micm::VectorMatrix; -template -using Group2VectorMatrix = micm::VectorMatrix; -template -using Group3VectorMatrix = micm::VectorMatrix; -template -using Group4VectorMatrix = micm::VectorMatrix; - -using Group1SparseVectorMatrix = micm::SparseMatrix>; -using Group2SparseVectorMatrix = micm::SparseMatrix>; -using Group3SparseVectorMatrix = micm::SparseMatrix>; -using Group4SparseVectorMatrix = micm::SparseMatrix>; - void run_solver(auto& solver) { auto state = solver.GetState(); @@ -145,17 +120,19 @@ void run_solver(auto& solver) { auto result = solver.Solve(time_step - elapsed_solve_time, state); elapsed_solve_time = result.final_time_; - state.variables_ = result.result_; } } } +template +using JitBuilder = micm::JitSolverBuilder; + TEST(JitRosenbrockSolver, AlphaMinusJacobian) { - testAlphaMinusJacobian<1, Group1VectorMatrix, Group1SparseVectorMatrix>(); - testAlphaMinusJacobian<2, Group2VectorMatrix, Group2SparseVectorMatrix>(); - testAlphaMinusJacobian<3, Group3VectorMatrix, Group3SparseVectorMatrix>(); - testAlphaMinusJacobian<4, Group4VectorMatrix, Group4SparseVectorMatrix>(); + testAlphaMinusJacobian(JitBuilder<1>(micm::JitRosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 1); + testAlphaMinusJacobian(JitBuilder<2>(micm::JitRosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 2); + testAlphaMinusJacobian(JitBuilder<3>(micm::JitRosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 3); + testAlphaMinusJacobian(JitBuilder<4>(micm::JitRosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 4); } TEST(JitRosenbrockSolver, MultipleInstances) @@ -169,26 +146,13 @@ TEST(JitRosenbrockSolver, MultipleInstances) auto chemical_system = solver_params.system_; auto reactions = solver_params.processes_; - auto solver_parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - - micm::JitRosenbrockSolver< - Group1VectorMatrix, - Group1SparseVectorMatrix, - micm::JitLinearSolver<1, Group1SparseVectorMatrix>, - micm::JitProcessSet<1>> - solver1(chemical_system, reactions, solver_parameters); - micm::JitRosenbrockSolver< - Group1VectorMatrix, - Group1SparseVectorMatrix, - micm::JitLinearSolver<1, Group1SparseVectorMatrix>, - micm::JitProcessSet<1>> - solver2(chemical_system, reactions, solver_parameters); - micm::JitRosenbrockSolver< - Group1VectorMatrix, - Group1SparseVectorMatrix, - micm::JitLinearSolver<1, Group1SparseVectorMatrix>, - micm::JitProcessSet<1>> - solver3(chemical_system, reactions, solver_parameters); + auto builder = JitBuilder<1>(micm::JitRosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(chemical_system) + .SetReactions(reactions); + + auto solver1 = builder.Build(); + auto solver2 = builder.Build(); + auto solver3 = builder.Build(); run_solver(solver1); run_solver(solver2); diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp index e795fc43c..a48355cde 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -1,5 +1,8 @@ #include #include +#include +#include +#include #include #include #include @@ -92,12 +95,12 @@ TEST(SolverBuilder, CanBuildRosenbrock) TEST(SolverBuilder, CanBuildJitRosenbrock) { - // auto jit_rosenbrock = micm::JitSolverBuilder() - // .SetSystem(the_system) - // .SetReactions(reactions) - // .SetNumberOfGridCells(1) - // .SolverParameters(micm::ThreeStageRosenbockSolverParameters{}) - // .Build(); + constexpr std::size_t L = 4; + auto jit_rosenbrock = micm::JitSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(L) + .Build(); } TEST(SolverBuilder, CanBuildCudaSolvers) diff --git a/test/unit/util/test_sparse_matrix_policy.hpp b/test/unit/util/test_sparse_matrix_policy.hpp index 10ec8c8b2..cfdf58642 100644 --- a/test/unit/util/test_sparse_matrix_policy.hpp +++ b/test/unit/util/test_sparse_matrix_policy.hpp @@ -162,27 +162,29 @@ MatrixPolicy testSingleBlockMatrix() .WithElement(3, 2) .WithElement(0, 1) .WithElement(2, 3) + .WithElement(1, 1) .WithElement(2, 1); // 0 X 0 0 - // 0 0 0 0 + // 0 X 0 0 // 0 X 0 X // 0 0 X 0 { auto row_ids = builder.RowIdsVector(); auto row_starts = builder.RowStartVector(); - EXPECT_EQ(builder.NumberOfElements(), 4); - EXPECT_EQ(row_ids.size(), 4); + EXPECT_EQ(builder.NumberOfElements(), 5); + EXPECT_EQ(row_ids.size(), 5); EXPECT_EQ(row_ids[0], 1); EXPECT_EQ(row_ids[1], 1); - EXPECT_EQ(row_ids[2], 3); - EXPECT_EQ(row_ids[3], 2); + EXPECT_EQ(row_ids[2], 1); + EXPECT_EQ(row_ids[3], 3); + EXPECT_EQ(row_ids[4], 2); EXPECT_EQ(row_starts.size(), 5); EXPECT_EQ(row_starts[0], 0); EXPECT_EQ(row_starts[1], 1); - EXPECT_EQ(row_starts[2], 1); - EXPECT_EQ(row_starts[3], 3); - EXPECT_EQ(row_starts[4], 4); + EXPECT_EQ(row_starts[2], 2); + EXPECT_EQ(row_starts[3], 4); + EXPECT_EQ(row_starts[4], 5); } MatrixPolicy matrix{ builder }; @@ -190,20 +192,25 @@ MatrixPolicy testSingleBlockMatrix() { auto& row_ids = matrix.RowIdsVector(); auto& row_starts = matrix.RowStartVector(); - EXPECT_EQ(row_ids.size(), 4); + EXPECT_EQ(row_ids.size(), 5); EXPECT_EQ(row_ids[0], 1); EXPECT_EQ(row_ids[1], 1); - EXPECT_EQ(row_ids[2], 3); - EXPECT_EQ(row_ids[3], 2); + EXPECT_EQ(row_ids[2], 1); + EXPECT_EQ(row_ids[3], 3); + EXPECT_EQ(row_ids[4], 2); EXPECT_EQ(row_starts.size(), 5); EXPECT_EQ(row_starts[0], 0); EXPECT_EQ(row_starts[1], 1); - EXPECT_EQ(row_starts[2], 1); - EXPECT_EQ(row_starts[3], 3); - EXPECT_EQ(row_starts[4], 4); + EXPECT_EQ(row_starts[2], 2); + EXPECT_EQ(row_starts[3], 4); + EXPECT_EQ(row_starts[4], 5); + + auto diagonal_ids = matrix.DiagonalIndices(0); + EXPECT_EQ(diagonal_ids.size(), 1); + EXPECT_EQ(diagonal_ids[0], matrix.VectorIndex(1, 1)); } - EXPECT_EQ(matrix.FlatBlockSize(), 4); + EXPECT_EQ(matrix.FlatBlockSize(), 5); EXPECT_EQ(matrix.IsZero(0, 1), false); EXPECT_EQ(matrix.IsZero(3, 2), false); @@ -235,7 +242,7 @@ MatrixPolicy testSingleBlockMatrix() }, std::system_error); EXPECT_THROW( - try { std::size_t elem = matrix.VectorIndex(1, 1); } catch (const std::system_error& e) { + try { std::size_t elem = matrix.VectorIndex(2, 2); } catch (const std::system_error& e) { EXPECT_EQ(e.code().value(), static_cast(MicmMatrixErrc::ZeroElementAccess)); throw; }, @@ -265,7 +272,7 @@ MatrixPolicy testSingleBlockMatrix() }, std::system_error); EXPECT_THROW( - try { matrix[0][1][1] = 2; } catch (const std::system_error& e) { + try { matrix[0][3][3] = 2; } catch (const std::system_error& e) { EXPECT_EQ(e.code().value(), static_cast(MicmMatrixErrc::ZeroElementAccess)); throw; }, @@ -281,9 +288,10 @@ MatrixPolicy testConstSingleBlockMatrix() .WithElement(3, 2) .WithElement(0, 1) .WithElement(2, 3) + .WithElement(1, 1) .WithElement(2, 1); // 0 X 0 0 - // 0 0 0 0 + // 0 X 0 0 // 0 X 0 X // 0 0 X 0 MatrixPolicy orig_matrix{ builder }; @@ -291,26 +299,31 @@ MatrixPolicy testConstSingleBlockMatrix() orig_matrix[0][3][2] = 42; orig_matrix[0][2][3] = 21; auto& row_ids = orig_matrix.RowIdsVector(); - EXPECT_EQ(row_ids.size(), 4); + EXPECT_EQ(row_ids.size(), 5); const MatrixPolicy matrix = orig_matrix; { auto& row_ids = matrix.RowIdsVector(); auto& row_starts = matrix.RowStartVector(); - EXPECT_EQ(row_ids.size(), 4); + EXPECT_EQ(row_ids.size(), 5); EXPECT_EQ(row_ids[0], 1); EXPECT_EQ(row_ids[1], 1); - EXPECT_EQ(row_ids[2], 3); - EXPECT_EQ(row_ids[3], 2); + EXPECT_EQ(row_ids[2], 1); + EXPECT_EQ(row_ids[3], 3); + EXPECT_EQ(row_ids[4], 2); EXPECT_EQ(row_starts.size(), 5); EXPECT_EQ(row_starts[0], 0); EXPECT_EQ(row_starts[1], 1); - EXPECT_EQ(row_starts[2], 1); - EXPECT_EQ(row_starts[3], 3); - EXPECT_EQ(row_starts[4], 4); + EXPECT_EQ(row_starts[2], 2); + EXPECT_EQ(row_starts[3], 4); + EXPECT_EQ(row_starts[4], 5); + + auto diagonal_ids = matrix.DiagonalIndices(0); + EXPECT_EQ(diagonal_ids.size(), 1); + EXPECT_EQ(diagonal_ids[0], matrix.VectorIndex(1, 1)); } - EXPECT_EQ(matrix.FlatBlockSize(), 4); + EXPECT_EQ(matrix.FlatBlockSize(), 5); EXPECT_EQ(matrix.IsZero(0, 1), false); EXPECT_EQ(matrix.IsZero(3, 2), false); @@ -340,7 +353,7 @@ MatrixPolicy testConstSingleBlockMatrix() }, std::system_error); EXPECT_THROW( - try { std::size_t elem = matrix.VectorIndex(1, 1); } catch (const std::system_error& e) { + try { std::size_t elem = matrix.VectorIndex(2, 2); } catch (const std::system_error& e) { EXPECT_EQ(e.code().value(), static_cast(MicmMatrixErrc::ZeroElementAccess)); throw; }, @@ -362,32 +375,34 @@ MatrixPolicy testMultiBlockMatrix() .WithElement(3, 2) .WithElement(0, 1) .WithElement(2, 3) + .WithElement(1, 1) .WithElement(2, 1) .InitialValue(24) .SetNumberOfBlocks(3); // 0 X 0 0 - // 0 0 0 0 + // 0 X 0 0 // 0 X 0 X // 0 0 X 0 auto row_ids = builder.RowIdsVector(); auto row_starts = builder.RowStartVector(); - EXPECT_EQ(builder.NumberOfElements(), 4 * 3); - EXPECT_EQ(row_ids.size(), 4); + EXPECT_EQ(builder.NumberOfElements(), 5 * 3); + EXPECT_EQ(row_ids.size(), 5); EXPECT_EQ(row_ids[0], 1); EXPECT_EQ(row_ids[1], 1); - EXPECT_EQ(row_ids[2], 3); - EXPECT_EQ(row_ids[3], 2); + EXPECT_EQ(row_ids[2], 1); + EXPECT_EQ(row_ids[3], 3); + EXPECT_EQ(row_ids[4], 2); EXPECT_EQ(row_starts.size(), 5); EXPECT_EQ(row_starts[0], 0); EXPECT_EQ(row_starts[1], 1); - EXPECT_EQ(row_starts[2], 1); - EXPECT_EQ(row_starts[3], 3); - EXPECT_EQ(row_starts[4], 4); + EXPECT_EQ(row_starts[2], 2); + EXPECT_EQ(row_starts[3], 4); + EXPECT_EQ(row_starts[4], 5); MatrixPolicy matrix{ builder }; - EXPECT_EQ(matrix.FlatBlockSize(), 4); + EXPECT_EQ(matrix.FlatBlockSize(), 5); EXPECT_EQ(matrix[0][2][1], 24); matrix[0][2][1] = 45; @@ -396,6 +411,16 @@ MatrixPolicy testMultiBlockMatrix() matrix[2][2][3] = 64; EXPECT_EQ(matrix[2][2][3], 64); + auto diagonal_ids = matrix.DiagonalIndices(0); + EXPECT_EQ(diagonal_ids.size(), 1); + EXPECT_EQ(diagonal_ids[0], matrix.VectorIndex(0, 1, 1)); + diagonal_ids = matrix.DiagonalIndices(1); + EXPECT_EQ(diagonal_ids.size(), 1); + EXPECT_EQ(diagonal_ids[0], matrix.VectorIndex(1, 1, 1)); + diagonal_ids = matrix.DiagonalIndices(2); + EXPECT_EQ(diagonal_ids.size(), 1); + EXPECT_EQ(diagonal_ids[0], matrix.VectorIndex(2, 1, 1)); + EXPECT_THROW( try { std::size_t elem = matrix.VectorIndex(0, 4, 2); } catch (const std::system_error& e) { EXPECT_EQ(e.code().value(), static_cast(MicmMatrixErrc::ElementOutOfRange)); @@ -415,7 +440,7 @@ MatrixPolicy testMultiBlockMatrix() }, std::system_error); EXPECT_THROW( - try { std::size_t elem = matrix.VectorIndex(1, 1, 1); } catch (const std::system_error& e) { + try { std::size_t elem = matrix.VectorIndex(1, 2, 2); } catch (const std::system_error& e) { EXPECT_EQ(e.code().value(), static_cast(MicmMatrixErrc::ZeroElementAccess)); throw; }, @@ -451,7 +476,7 @@ MatrixPolicy testMultiBlockMatrix() }, std::system_error); EXPECT_THROW( - try { matrix[0][1][1] = 2; } catch (const std::system_error& e) { + try { matrix[0][3][3] = 2; } catch (const std::system_error& e) { EXPECT_EQ(e.code().value(), static_cast(MicmMatrixErrc::ZeroElementAccess)); throw; }, diff --git a/test/unit/util/test_sparse_matrix_standard_ordering.cpp b/test/unit/util/test_sparse_matrix_standard_ordering.cpp index 805e93a0f..315780d98 100644 --- a/test/unit/util/test_sparse_matrix_standard_ordering.cpp +++ b/test/unit/util/test_sparse_matrix_standard_ordering.cpp @@ -28,15 +28,15 @@ TEST(SparseMatrix, SingleBlockMatrix) { std::size_t elem = matrix.VectorIndex(3, 2); - EXPECT_EQ(elem, 3); + EXPECT_EQ(elem, 4); matrix.AsVector()[elem] = 42; - EXPECT_EQ(matrix.AsVector()[3], 42); + EXPECT_EQ(matrix.AsVector()[4], 42); } { std::size_t elem = matrix.VectorIndex(2, 3); - EXPECT_EQ(elem, 2); + EXPECT_EQ(elem, 3); matrix.AsVector()[elem] = 21; - EXPECT_EQ(matrix.AsVector()[2], 21); + EXPECT_EQ(matrix.AsVector()[3], 21); } } @@ -45,13 +45,13 @@ TEST(SparseMatrix, ConstSingleBlockMatrix) auto matrix = testConstSingleBlockMatrix(); { std::size_t elem = matrix.VectorIndex(3, 2); - EXPECT_EQ(elem, 3); - EXPECT_EQ(matrix.AsVector()[3], 42); + EXPECT_EQ(elem, 4); + EXPECT_EQ(matrix.AsVector()[4], 42); } { std::size_t elem = matrix.VectorIndex(2, 3); - EXPECT_EQ(elem, 2); - EXPECT_EQ(matrix.AsVector()[2], 21); + EXPECT_EQ(elem, 3); + EXPECT_EQ(matrix.AsVector()[3], 21); } } @@ -61,15 +61,15 @@ TEST(SparseMatrix, MultiBlockMatrix) { std::size_t elem = matrix.VectorIndex(0, 2, 3); - EXPECT_EQ(elem, 2); + EXPECT_EQ(elem, 3); matrix.AsVector()[elem] = 21; - EXPECT_EQ(matrix.AsVector()[2], 21); + EXPECT_EQ(matrix.AsVector()[3], 21); } { std::size_t elem = matrix.VectorIndex(2, 2, 1); - EXPECT_EQ(elem, 9); + EXPECT_EQ(elem, 12); matrix.AsVector()[elem] = 31; - EXPECT_EQ(matrix.AsVector()[9], 31); + EXPECT_EQ(matrix.AsVector()[12], 31); } } diff --git a/test/unit/util/test_sparse_matrix_vector_ordering.cpp b/test/unit/util/test_sparse_matrix_vector_ordering.cpp index 4d351af8c..68be6fec1 100644 --- a/test/unit/util/test_sparse_matrix_vector_ordering.cpp +++ b/test/unit/util/test_sparse_matrix_vector_ordering.cpp @@ -26,18 +26,18 @@ TEST(SparseVectorMatrix, SingleBlockMatrix) { std::size_t elem = matrix.VectorIndex(3, 2); - EXPECT_EQ(elem, 12); + EXPECT_EQ(elem, 16); matrix.AsVector()[elem] = 42; - EXPECT_EQ(matrix.AsVector()[12], 42); + EXPECT_EQ(matrix.AsVector()[16], 42); } { std::size_t elem = matrix.VectorIndex(2, 3); - EXPECT_EQ(elem, 8); + EXPECT_EQ(elem, 12); matrix.AsVector()[elem] = 21; - EXPECT_EQ(matrix.AsVector()[8], 21); + EXPECT_EQ(matrix.AsVector()[12], 21); } EXPECT_EQ(matrix.GroupVectorSize(), 4); - EXPECT_EQ(matrix.GroupSize(matrix.FlatBlockSize()), 4 * 4); + EXPECT_EQ(matrix.GroupSize(matrix.FlatBlockSize()), 5 * 4); EXPECT_EQ(matrix.NumberOfGroups(1), 1); } @@ -47,16 +47,16 @@ TEST(SparseVectorMatrix, ConstSingleBlockMatrix) { std::size_t elem = matrix.VectorIndex(3, 2); - EXPECT_EQ(elem, 6); - EXPECT_EQ(matrix.AsVector()[6], 42); + EXPECT_EQ(elem, 8); + EXPECT_EQ(matrix.AsVector()[8], 42); } { std::size_t elem = matrix.VectorIndex(2, 3); - EXPECT_EQ(elem, 4); - EXPECT_EQ(matrix.AsVector()[4], 21); + EXPECT_EQ(elem, 6); + EXPECT_EQ(matrix.AsVector()[6], 21); } EXPECT_EQ(matrix.GroupVectorSize(), 2); - EXPECT_EQ(matrix.GroupSize(matrix.FlatBlockSize()), 2 * 4); + EXPECT_EQ(matrix.GroupSize(matrix.FlatBlockSize()), 2 * 5); EXPECT_EQ(matrix.NumberOfGroups(1), 1); } @@ -66,17 +66,17 @@ TEST(SparseVectorMatrix, MultiBlockMatrix) { std::size_t elem = matrix.VectorIndex(0, 2, 3); - EXPECT_EQ(elem, 4); + EXPECT_EQ(elem, 6); matrix.AsVector()[elem] = 21; - EXPECT_EQ(matrix.AsVector()[4], 21); + EXPECT_EQ(matrix.AsVector()[6], 21); } { std::size_t elem = matrix.VectorIndex(2, 2, 1); - EXPECT_EQ(elem, 10); + EXPECT_EQ(elem, 14); matrix.AsVector()[elem] = 31; - EXPECT_EQ(matrix.AsVector()[10], 31); + EXPECT_EQ(matrix.AsVector()[14], 31); } EXPECT_EQ(matrix.GroupVectorSize(), 2); - EXPECT_EQ(matrix.GroupSize(matrix.FlatBlockSize()), 2 * 4); + EXPECT_EQ(matrix.GroupSize(matrix.FlatBlockSize()), 2 * 5); EXPECT_EQ(matrix.NumberOfGroups(4), 2); }