diff --git a/docker/Dockerfile.intel b/docker/Dockerfile.intel index cb41cef7c..70fb8e9b8 100644 --- a/docker/Dockerfile.intel +++ b/docker/Dockerfile.intel @@ -7,7 +7,6 @@ RUN apt update \ ca-certificates \ cmake \ cmake-curses-gui \ - curl \ libcurl4-openssl-dev \ libhdf5-dev \ m4 \ diff --git a/docker/Dockerfile.llvm b/docker/Dockerfile.llvm index 33b3272fa..005f24b1d 100644 --- a/docker/Dockerfile.llvm +++ b/docker/Dockerfile.llvm @@ -10,6 +10,7 @@ RUN dnf -y update \ make \ zlib-devel \ llvm-devel \ + openmpi-devel \ valgrind \ && dnf clean all @@ -24,6 +25,7 @@ RUN mkdir /build \ -D MICM_ENABLE_CLANG_TIDY:BOOL=FALSE \ -D MICM_ENABLE_LLVM:BOOL=TRUE \ -D MICM_ENABLE_MEMCHECK:BOOL=TRUE \ + -D MICM_ENABLE_OPENMP:BOOL=TRUE \ ../micm \ && make install -j 8 diff --git a/examples/profile_example.cpp b/examples/profile_example.cpp index 48f82d21c..25f380e5b 100644 --- a/examples/profile_example.cpp +++ b/examples/profile_example.cpp @@ -29,7 +29,7 @@ namespace fs = std::filesystem; using namespace micm; -template class MatrixType, template class SparseMatrixType> +template class MatrixType, class SparseMatrixType> int Run(const char* filepath, const char* initial_conditions, const std::string& matrix_ordering_type) { using SolverType = RosenbrockSolver; @@ -114,24 +114,9 @@ int Run(const char* filepath, const char* initial_conditions, const std::string& return 0; } -template -using SparseMatrixParam = micm::SparseMatrix; -template -using Vector1MatrixParam = micm::VectorMatrix; -template -using Vector10MatrixParam = micm::VectorMatrix; -template -using Vector100MatrixParam = micm::VectorMatrix; -template +template using Vector1000MatrixParam = micm::VectorMatrix; -template -using Vector1SparseMatrixParam = micm::SparseMatrix>; -template -using Vector10SparseMatrixParam = micm::SparseMatrix>; -template -using Vector100SparseMatrixParam = micm::SparseMatrix>; -template -using Vector1000SparseMatrixParam = micm::SparseMatrix>; +using Vector1000SparseMatrixParam = micm::SparseMatrix>; int main(const int argc, const char* argv[]) { diff --git a/include/micm/jit/jit_compiler.hpp b/include/micm/jit/jit_compiler.hpp index 66e0c7cba..7637dc3b9 100644 --- a/include/micm/jit/jit_compiler.hpp +++ b/include/micm/jit/jit_compiler.hpp @@ -40,7 +40,8 @@ enum class MicmJitErrc { InvalidMatrix = MICM_JIT_ERROR_CODE_INVALID_MATRIX, - MissingJitFunction = MICM_JIT_ERROR_CODE_MISSING_JIT_FUNCTION + MissingJitFunction = MICM_JIT_ERROR_CODE_MISSING_JIT_FUNCTION, + FailedToBuild = MICM_JIT_ERROR_CODE_FAILED_TO_BUILD }; namespace std @@ -84,6 +85,7 @@ inline std::error_code make_error_code(MicmJitErrc e) namespace micm { + // a singleton class class JitCompiler { private: @@ -99,23 +101,23 @@ namespace micm llvm::orc::JITDylib &main_lib_; public: - JitCompiler( - std::unique_ptr execution_session, - llvm::orc::JITTargetMachineBuilder machine_builder, - llvm::DataLayout data_layout) - : execution_session_(std::move(execution_session)), - data_layout_(std::move(data_layout)), - mangle_(*this->execution_session_, this->data_layout_), - object_layer_(*this->execution_session_, []() { return std::make_unique(); }), - compile_layer_( - *this->execution_session_, - object_layer_, - std::make_unique(std::move(machine_builder))), - optimize_layer_(*this->execution_session_, compile_layer_, OptimizeModule), - main_lib_(this->execution_session_->createBareJITDylib("
")) + // Delete the copy constructor and assignment operator + JitCompiler(const JitCompiler &) = delete; + JitCompiler &operator=(const JitCompiler &) = delete; + + static JitCompiler &GetInstance() { - main_lib_.addGenerator( - llvm::cantFail(llvm::orc::DynamicLibrarySearchGenerator::GetForCurrentProcess(data_layout_.getGlobalPrefix()))); + static std::unique_ptr instance; + if (!instance) + { + auto expectedInstance = Create(); + if (!expectedInstance) + { + throw std::system_error(make_error_code(MicmJitErrc::FailedToBuild)); + } + instance = std::move(*expectedInstance); + } + return *instance; } ~JitCompiler() @@ -124,28 +126,6 @@ namespace micm execution_session_->reportError(std::move(Err)); } - static llvm::Expected> Create() - { - llvm::InitializeNativeTarget(); - llvm::InitializeNativeTargetAsmPrinter(); - llvm::InitializeNativeTargetAsmParser(); - - auto EPC = llvm::orc::SelfExecutorProcessControl::Create(); - if (!EPC) - return EPC.takeError(); - - auto execution_session = std::make_unique(std::move(*EPC)); - - llvm::orc::JITTargetMachineBuilder machine_builder(execution_session->getExecutorProcessControl().getTargetTriple()); - - auto data_layout = machine_builder.getDefaultDataLayoutForTarget(); - if (!data_layout) - return data_layout.takeError(); - - return std::make_shared( - std::move(execution_session), std::move(machine_builder), std::move(*data_layout)); - } - const llvm::DataLayout &GetDataLayout() const { return data_layout_; @@ -171,6 +151,47 @@ namespace micm } private: + static llvm::Expected> Create() + { + llvm::InitializeNativeTarget(); + llvm::InitializeNativeTargetAsmPrinter(); + llvm::InitializeNativeTargetAsmParser(); + + auto EPC = llvm::orc::SelfExecutorProcessControl::Create(); + if (!EPC) + return EPC.takeError(); + + auto execution_session = std::make_unique(std::move(*EPC)); + + llvm::orc::JITTargetMachineBuilder machine_builder(execution_session->getExecutorProcessControl().getTargetTriple()); + + auto data_layout = machine_builder.getDefaultDataLayoutForTarget(); + if (!data_layout) + return data_layout.takeError(); + + return llvm::Expected>(std::unique_ptr( + new JitCompiler(std::move(execution_session), std::move(machine_builder), std::move(*data_layout)))); + } + + JitCompiler( + std::unique_ptr execution_session, + llvm::orc::JITTargetMachineBuilder machine_builder, + llvm::DataLayout data_layout) + : execution_session_(std::move(execution_session)), + data_layout_(std::move(data_layout)), + mangle_(*this->execution_session_, this->data_layout_), + object_layer_(*this->execution_session_, []() { return std::make_unique(); }), + compile_layer_( + *this->execution_session_, + object_layer_, + std::make_unique(std::move(machine_builder))), + optimize_layer_(*this->execution_session_, compile_layer_, OptimizeModule), + main_lib_(this->execution_session_->createBareJITDylib("
")) + { + main_lib_.addGenerator( + llvm::cantFail(llvm::orc::DynamicLibrarySearchGenerator::GetForCurrentProcess(data_layout_.getGlobalPrefix()))); + } + static llvm::Expected OptimizeModule( llvm::orc::ThreadSafeModule threadsafe_module, const llvm::orc::MaterializationResponsibility &responsibility) diff --git a/include/micm/jit/jit_function.hpp b/include/micm/jit/jit_function.hpp index e14273f16..2b970e45a 100644 --- a/include/micm/jit/jit_function.hpp +++ b/include/micm/jit/jit_function.hpp @@ -77,7 +77,7 @@ namespace micm { bool generated_ = false; std::string name_; - std::shared_ptr compiler_; + JitCompiler* compiler_; public: std::unique_ptr context_; @@ -91,7 +91,7 @@ namespace micm JitFunction() = delete; friend class JitFunctionBuilder; - static JitFunctionBuilder Create(std::shared_ptr compiler); + static JitFunctionBuilder Create(); JitFunction(JitFunctionBuilder& function_builder); /// @brief Generates the function @@ -138,7 +138,7 @@ namespace micm class JitFunctionBuilder { - std::shared_ptr compiler_; + JitCompiler* compiler_; std::string name_; std::vector> arguments_; JitType return_type_{ JitType::Void }; @@ -146,15 +146,15 @@ namespace micm public: JitFunctionBuilder() = delete; - JitFunctionBuilder(std::shared_ptr compiler); + JitFunctionBuilder(JitCompiler& compiler); JitFunctionBuilder& SetName(const std::string& name); JitFunctionBuilder& SetArguments(const std::vector>& arguments); JitFunctionBuilder& SetReturnType(JitType type); }; - inline JitFunctionBuilder JitFunction::Create(std::shared_ptr compiler) + inline JitFunctionBuilder JitFunction::Create() { - return JitFunctionBuilder{ compiler }; + return JitFunctionBuilder{ JitCompiler::GetInstance() }; } JitFunction::JitFunction(JitFunctionBuilder& function_builder) @@ -296,8 +296,8 @@ namespace micm return TmpB.CreateAlloca(type, 0, var_name.c_str()); } - inline JitFunctionBuilder::JitFunctionBuilder(std::shared_ptr compiler) - : compiler_(compiler){}; + inline JitFunctionBuilder::JitFunctionBuilder(JitCompiler& compiler) + : compiler_(&compiler){}; inline JitFunctionBuilder& JitFunctionBuilder::SetName(const std::string& name) { diff --git a/include/micm/process/cuda_process_set.hpp b/include/micm/process/cuda_process_set.hpp index 462eca7e2..f678ac30d 100644 --- a/include/micm/process/cuda_process_set.hpp +++ b/include/micm/process/cuda_process_set.hpp @@ -33,17 +33,17 @@ namespace micm template void SetJacobianFlatIds(const SparseMatrix& matrix); - template typename MatrixPolicy> - requires(CudaMatrix>&& VectorizableDense>) void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; - - template typename MatrixPolicy> - requires(!CudaMatrix>) void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; + template + requires(CudaMatrix&& VectorizableDense) void AddForcingTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + MatrixPolicy& forcing) const; + + template + requires(!CudaMatrix) void AddForcingTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + MatrixPolicy& forcing) const; template requires( @@ -101,12 +101,12 @@ namespace micm micm::cuda::CopyJacobiFlatId(hoststruct, this->devstruct_); } - template class MatrixPolicy> - requires(CudaMatrix>&& VectorizableDense>) inline void CudaProcessSet:: + template + requires(CudaMatrix&& VectorizableDense) inline void CudaProcessSet:: AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + MatrixPolicy& forcing) const { auto forcing_param = forcing.AsDeviceParam(); // we need to update forcing so it can't be constant and must be an lvalue micm::cuda::AddForcingTermsKernelDriver( @@ -114,11 +114,11 @@ namespace micm } // call the function from the base class - template class MatrixPolicy> - requires(!CudaMatrix>) inline void CudaProcessSet::AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const + template + requires(!CudaMatrix) inline void CudaProcessSet::AddForcingTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + MatrixPolicy& forcing) const { AddForcingTerms(rate_constants, state_variables, forcing); } diff --git a/include/micm/process/jit_process_set.hpp b/include/micm/process/jit_process_set.hpp index 0a92e9780..aaeabbd13 100644 --- a/include/micm/process/jit_process_set.hpp +++ b/include/micm/process/jit_process_set.hpp @@ -4,7 +4,6 @@ */ #pragma once -#include #include #include #include @@ -19,7 +18,6 @@ namespace micm template class JitProcessSet : public ProcessSet { - std::shared_ptr compiler_; llvm::orc::ResourceTrackerSP forcing_function_resource_tracker_; void (*forcing_function_)(const double *, const double *, double *); llvm::orc::ResourceTrackerSP jacobian_function_resource_tracker_; @@ -34,11 +32,9 @@ namespace micm JitProcessSet() = default; /// @brief Create a JITed process set calculator for a given set of processes - /// @param compiler JIT compiler /// @param processes Processes to create calculator for /// @param variable_map A mapping of species names to concentration index JitProcessSet( - std::shared_ptr compiler, const std::vector &processes, const std::map &variable_map); @@ -53,11 +49,11 @@ namespace micm /// @param rate_constants Current values for the process rate constants (grid cell, process) /// @param state_variables Current state variable values (grid cell, state variable) /// @param forcing Forcing terms for each state variable (grid cell, state variable) - template class MatrixPolicy> + template void AddForcingTerms( - const MatrixPolicy &rate_constants, - const MatrixPolicy &state_variables, - MatrixPolicy &forcing) const; + const MatrixPolicy &rate_constants, + const MatrixPolicy &state_variables, + MatrixPolicy &forcing) const; /// @brief Subtracts Jacobian terms for the set of processes for the current conditions /// @param rate_constants Current values for the process rate constants (grid cell, process) @@ -81,7 +77,6 @@ namespace micm template inline JitProcessSet::JitProcessSet(JitProcessSet &&other) : ProcessSet(std::move(other)), - compiler_(std::move(other.compiler_)), forcing_function_resource_tracker_(std::move(other.forcing_function_resource_tracker_)), forcing_function_(std::move(other.forcing_function_)), jacobian_function_resource_tracker_(std::move(other.jacobian_function_resource_tracker_)), @@ -95,7 +90,6 @@ namespace micm inline JitProcessSet &JitProcessSet::operator=(JitProcessSet &&other) { ProcessSet::operator=(std::move(other)); - compiler_ = std::move(other.compiler_); forcing_function_resource_tracker_ = std::move(other.forcing_function_resource_tracker_); forcing_function_ = std::move(other.forcing_function_); jacobian_function_resource_tracker_ = std::move(other.jacobian_function_resource_tracker_); @@ -107,11 +101,9 @@ namespace micm template inline JitProcessSet::JitProcessSet( - std::shared_ptr compiler, const std::vector &processes, const std::map &variable_map) - : ProcessSet(processes, variable_map), - compiler_(compiler) + : ProcessSet(processes, variable_map) { forcing_function_ = NULL; jacobian_function_ = NULL; @@ -122,7 +114,7 @@ namespace micm void JitProcessSet::GenerateForcingFunction() { std::string function_name = "add_forcing_terms_" + GenerateRandomString(); - JitFunction func = JitFunction::Create(compiler_) + JitFunction func = JitFunction::Create() .SetName(function_name) .SetArguments({ { "rate constants", JitType::DoublePtr }, { "state variables", JitType::DoublePtr }, @@ -224,7 +216,7 @@ namespace micm void JitProcessSet::GenerateJacobianFunction(const SparseMatrix> &matrix) { std::string function_name = "subtract_jacobian_terms_" + GenerateRandomString(); - JitFunction func = JitFunction::Create(compiler_) + JitFunction func = JitFunction::Create() .SetName(function_name) .SetArguments({ { "rate constants", JitType::DoublePtr }, { "state variables", JitType::DoublePtr }, @@ -335,11 +327,11 @@ namespace micm } template - template class MatrixPolicy> + template void JitProcessSet::AddForcingTerms( - const MatrixPolicy &rate_constants, - const MatrixPolicy &state_variables, - MatrixPolicy &forcing) const + const MatrixPolicy &rate_constants, + const MatrixPolicy &state_variables, + MatrixPolicy &forcing) const { forcing_function_(rate_constants.AsVector().data(), state_variables.AsVector().data(), forcing.AsVector().data()); } diff --git a/include/micm/process/process.hpp b/include/micm/process/process.hpp index 165f9357c..f190315e2 100644 --- a/include/micm/process/process.hpp +++ b/include/micm/process/process.hpp @@ -85,14 +85,14 @@ namespace micm /// @brief Update the solver state rate constants /// @param processes The set of processes being solved /// @param state The solver state to update - template class MatrixPolicy, template class SparseMatrixPolicy> - requires(!VectorizableDense>) static void UpdateState( + template + requires(!VectorizableDense) static void UpdateState( const std::vector& processes, - State& state); - template class MatrixPolicy, template class SparseMatrixPolicy> - requires(VectorizableDense>) static void UpdateState( + State& state); + template + requires(VectorizableDense) static void UpdateState( const std::vector& processes, - State& state); + State& state); friend class ProcessBuilder; static ProcessBuilder Create(); @@ -148,10 +148,10 @@ namespace micm ProcessBuilder& SetPhase(const Phase& phase); }; - template class MatrixPolicy, template class SparseMatrixPolicy> - requires(!VectorizableDense>) void Process::UpdateState( + template + requires(!VectorizableDense) void Process::UpdateState( const std::vector& processes, - State& state) + State& state) { MICM_PROFILE_FUNCTION(); @@ -173,10 +173,10 @@ namespace micm } } - template class MatrixPolicy, template class SparseMatrixPolicy> - requires(VectorizableDense>) void Process::UpdateState( + template + requires(VectorizableDense) void Process::UpdateState( const std::vector& processes, - State& state) + State& state) { MICM_PROFILE_FUNCTION(); const auto& v_custom_parameters = state.custom_rate_parameters_.AsVector(); diff --git a/include/micm/process/process_set.hpp b/include/micm/process/process_set.hpp index dca9eb6ea..f52424aa6 100644 --- a/include/micm/process/process_set.hpp +++ b/include/micm/process/process_set.hpp @@ -91,33 +91,32 @@ namespace micm /// @param rate_constants Current values for the process rate constants (grid cell, process) /// @param state_variables Current state variable values (grid cell, state variable) /// @param forcing Forcing terms for each state variable (grid cell, state variable) - template typename MatrixPolicy> - requires(!VectorizableDense>) void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; - template typename MatrixPolicy> - requires VectorizableDense> + template + requires(!VectorizableDense) void AddForcingTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const; + template + requires VectorizableDense void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const; /// @brief Subtract Jacobian terms for the set of processes for the current conditions /// @param rate_constants Current values for the process rate constants (grid cell, process) /// @param state_variables Current state variable values (grid cell, state variable) /// @param jacobian Jacobian matrix for the system (grid cell, dependent variable, independent variable) - // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(!VectorizableDense || !VectorizableSparse) void SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + template + requires(!VectorizableDense || !VectorizableSparse) void SubtractJacobianTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(VectorizableDense&& VectorizableSparse) void SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + template + requires(VectorizableDense&& VectorizableSparse) void SubtractJacobianTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; /// @brief Returns the set of species used in a set of processes @@ -217,11 +216,11 @@ namespace micm } } - template typename MatrixPolicy> - requires(!VectorizableDense>) inline void ProcessSet::AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const + template + requires(!VectorizableDense) inline void ProcessSet::AddForcingTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const { MICM_PROFILE_FUNCTION(); @@ -261,12 +260,12 @@ namespace micm } }; - template typename MatrixPolicy> - requires VectorizableDense> + template + requires VectorizableDense inline void ProcessSet::AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const { MICM_PROFILE_FUNCTION(); @@ -306,12 +305,11 @@ namespace micm } // Forming the Jacobian matrix "J" and returning "-J" to be consistent with the CUDA implementation - // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(!VectorizableDense || !VectorizableSparse) inline void ProcessSet:: + template + requires(!VectorizableDense || !VectorizableSparse) inline void ProcessSet:: SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const { MICM_PROFILE_FUNCTION(); @@ -362,12 +360,11 @@ namespace micm } // Forming the Jacobian matrix "J" and returning "-J" to be consistent with the CUDA implementation - // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(VectorizableDense&& VectorizableSparse) inline void ProcessSet:: + template + requires(VectorizableDense&& VectorizableSparse) inline void ProcessSet:: SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const { MICM_PROFILE_FUNCTION(); diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index fb53e4f54..8288a51bd 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/include/micm/solver/backward_euler.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -28,34 +29,40 @@ namespace micm { /// @brief An implementation of the fully implicit backward euler method + template class BackwardEuler { + BackwardEulerSolverParameters parameters_; + LinearSolverPolicy linear_solver_; + ProcessSetPolicy process_set_; + std::vector jacobian_diagonal_elements_; + std::vector processes_; + public: /// @brief Default constructor - BackwardEuler(); - - /// @brief Builds a backward eulder solver for the given system and processes - /// @param system The chemical system to create the solver for - /// @param processes The collection of chemical processes that will be applied during solving - BackwardEuler(const System& system, const std::vector& processes); + BackwardEuler( + BackwardEulerSolverParameters parameters, + LinearSolverPolicy linear_solver, + ProcessSetPolicy process_set, + std::vector jacobian_diagonal_elements, + std::vector& processes + ) : parameters_(parameters), + linear_solver_(linear_solver), + process_set_(process_set), + jacobian_diagonal_elements_(jacobian_diagonal_elements), + processes_(processes) + { + } virtual ~BackwardEuler() = default; /// @brief Advances the given step over the specified time step /// @param time_step Time [s] to advance the state by /// @param state The state to advance - /// @param linear_solver The linear solver to use - /// @param process_set The process set to use - /// @param processes The collection of chemical processes that will be applied during solving - /// @param jacobian_diagonal_elements The diagonal elements of the jacobian matrix - /// @return A struct containing results and a status code + /// @return Nothing, but the state is updated void Solve( double time_step, - auto& state, - auto linear_solver, - auto process_set, - const std::vector& processes, - auto jacobian_diagonal_elements); + auto& state); }; } // namespace micm diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index 9efd98d81..7bae79d9f 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -44,17 +44,10 @@ inline std::error_code make_error_code(MicmBackwardEulerErrc e) namespace micm { - inline BackwardEuler::BackwardEuler() - { - } - - inline void BackwardEuler::Solve( + template + inline void BackwardEuler::Solve( double time_step, - auto& state, - auto linear_solver, - auto process_set, - const std::vector& processes, - auto jacobian_diagonal_elements) + auto& state) { // A fully implicit euler implementation is given by the following equation: // y_{n+1} = y_n + H * f(t_{n+1}, y_{n+1}) @@ -67,21 +60,24 @@ namespace micm // if that fails, try H = H/10 once // if that fails, accept the current H but do not update the Yn vector + + double tolerance = parameters_.absolute_tolerance_[0]; + double small = parameters_.small; + std::size_t max_iter = parameters_.max_number_of_steps_; + const auto time_step_reductions = parameters_.time_step_reductions; + double H = time_step; double t = 0.0; - double tol = 1e-4; - double small = 1.0e-40; - bool singular = false; - std::size_t max_iter = 11; std::size_t n_successful_integrations = 0; std::size_t n_convergence_failures = 0; - const std::array time_step_reductions{ 0.5, 0.5, 0.5, 0.5, 0.1 }; + + bool singular = false; auto Yn = state.variables_; auto Yn1 = state.variables_; auto forcing = state.variables_; - Process::UpdateState(processes, state); + Process::UpdateState(processes_, state); while (t < time_step) { @@ -97,11 +93,11 @@ namespace micm // so we can use Yn1 to calculate the forcing and jacobian // calculate forcing std::fill(forcing.AsVector().begin(), forcing.AsVector().end(), 0.0); - process_set.AddForcingTerms(state.rate_constants_, Yn1, forcing); + process_set_.AddForcingTerms(state.rate_constants_, Yn1, forcing); // calculate jacobian std::fill(state.jacobian_.AsVector().begin(), state.jacobian_.AsVector().end(), 0.0); - process_set.SubtractJacobianTerms(state.rate_constants_, Yn1, state.jacobian_); + process_set_.SubtractJacobianTerms(state.rate_constants_, Yn1, state.jacobian_); // subtract the inverse of the time step from the diagonal // TODO: handle vectorized jacobian matrix @@ -112,7 +108,7 @@ namespace micm for (std::size_t i_block = 0; i_block < state.jacobian_.NumberOfBlocks(); ++i_block) { auto jacobian_vector = std::next(state.jacobian_.AsVector().begin(), i_block * state.jacobian_.FlatBlockSize()); - for (const auto& i_elem : jacobian_diagonal_elements) + for (const auto& i_elem : jacobian_diagonal_elements_) jacobian_vector[i_elem] -= 1 / H; } @@ -120,7 +116,7 @@ namespace micm // (y_{n+1} - y_n) / H = f(t_{n+1}, y_{n+1}) // try to find the root by factoring and solving the linear system - linear_solver.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular); + linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular); auto yn1_iter = Yn1.begin(); auto yn_iter = Yn.begin(); @@ -135,7 +131,7 @@ namespace micm // the result of the linear solver will be stored in forcing // this represents the change in the solution - linear_solver.Solve(forcing, forcing, state.lower_matrix_, state.upper_matrix_); + linear_solver_.Solve(forcing, forcing, state.lower_matrix_, state.upper_matrix_); // solution_blk in camchem // Yn1 = Yn1 + residual; @@ -160,7 +156,7 @@ namespace micm do { // changes that are much smaller than the tolerance are negligible and we assume can be accepted - converged = (std::abs(*forcing_iter) <= small) || (std::abs(*forcing_iter) <= tol * std::abs(*yn1_iter)); + converged = (std::abs(*forcing_iter) <= small) || (std::abs(*forcing_iter) <= tolerance * std::abs(*yn1_iter)); ++forcing_iter, ++yn1_iter; } while (converged && forcing_iter != forcing.end()); diff --git a/include/micm/solver/backward_euler_solver_parameters.hpp b/include/micm/solver/backward_euler_solver_parameters.hpp new file mode 100644 index 000000000..5aa211c0a --- /dev/null +++ b/include/micm/solver/backward_euler_solver_parameters.hpp @@ -0,0 +1,24 @@ +/* Copyright (C) 2023-2024 National Center for Atmospheric Research + * + * SPDX-License-Identifier: Apache-2.0 + */ +#pragma once + +#include +#include +#include +#include + +namespace micm +{ + + /// @brief BackwardEuler solver parameters + struct BackwardEulerSolverParameters + { + std::vector absolute_tolerance_; + double small { 1.0e-40 }; + size_t max_number_of_steps_{ 11 }; + std::array time_step_reductions{ 0.5, 0.5, 0.5, 0.5, 0.1 }; + }; + +} // namespace micm \ No newline at end of file diff --git a/include/micm/solver/cuda_linear_solver.hpp b/include/micm/solver/cuda_linear_solver.hpp index a5df8e81a..bba93ab83 100644 --- a/include/micm/solver/cuda_linear_solver.hpp +++ b/include/micm/solver/cuda_linear_solver.hpp @@ -15,8 +15,8 @@ namespace micm { - template class SparseMatrixPolicy, class LuDecompositionPolicy = CudaLuDecomposition> - class CudaLinearSolver : public LinearSolver + template + class CudaLinearSolver : public LinearSolver { public: /// This is an instance of struct "LinearSolverParam" that holds @@ -32,17 +32,17 @@ namespace micm /// in this case, we will use the CudaLuDecomposition specified at line 15; /// See line 62 of "linear_solver.inl" for more details about how /// this lamda function works; - CudaLinearSolver(const SparseMatrixPolicy& matrix, T initial_value) - : CudaLinearSolver( + CudaLinearSolver(const SparseMatrixPolicy& matrix, typename SparseMatrixPolicy::value_type initial_value) + : CudaLinearSolver( matrix, initial_value, - [&](const SparseMatrixPolicy& m) -> LuDecompositionPolicy { return LuDecompositionPolicy(m); }){}; + [&](const SparseMatrixPolicy& m) -> LuDecompositionPolicy { return LuDecompositionPolicy(m); }){}; CudaLinearSolver( - const SparseMatrixPolicy& matrix, - T initial_value, - const std::function&)> create_lu_decomp) - : LinearSolver(matrix, initial_value, create_lu_decomp) + const SparseMatrixPolicy& matrix, + typename SparseMatrixPolicy::value_type initial_value, + const std::function create_lu_decomp) + : LinearSolver(matrix, initial_value, create_lu_decomp) { LinearSolverParam hoststruct; @@ -68,25 +68,24 @@ namespace micm micm::cuda::FreeConstData(this->devstruct_); }; - template class MatrixPolicy> + template requires( - CudaMatrix>&& CudaMatrix>&& VectorizableDense>&& - VectorizableSparse>) void Solve(const MatrixPolicy& b, MatrixPolicy& x, const SparseMatrixPolicy& L, const SparseMatrixPolicy& U) + CudaMatrix&& CudaMatrix&& VectorizableDense&& + VectorizableSparse) void Solve(const MatrixPolicy& b, MatrixPolicy& x, const SparseMatrixPolicy& L, const SparseMatrixPolicy& U) const { auto x_param = x.AsDeviceParam(); // we need to update x so it can't be constant and must be an lvalue micm::cuda::SolveKernelDriver(b.AsDeviceParam(), x_param, L.AsDeviceParam(), U.AsDeviceParam(), this->devstruct_); }; - template class MatrixPolicy> - requires(!CudaMatrix> && !CudaMatrix>) void Solve( - const MatrixPolicy& b, - MatrixPolicy& x, - const SparseMatrixPolicy& L, - const SparseMatrixPolicy& U) const + template + requires(!CudaMatrix && !CudaMatrix) void Solve( + const MatrixPolicy& b, + MatrixPolicy& x, + const SparseMatrixPolicy& L, + const SparseMatrixPolicy& U) const { - LinearSolver::template Solve(b, x, L, U); + LinearSolver::template Solve(b, x, L, U); }; }; } // namespace micm diff --git a/include/micm/solver/cuda_lu_decomposition.hpp b/include/micm/solver/cuda_lu_decomposition.hpp index ea0ccee92..38e659f92 100644 --- a/include/micm/solver/cuda_lu_decomposition.hpp +++ b/include/micm/solver/cuda_lu_decomposition.hpp @@ -28,9 +28,9 @@ namespace micm /// This is the overloaded constructor that takes one argument called "matrix"; /// We need to specify the type (e.g., double, int, etc) and /// ordering (e.g., vector-stored, non-vector-stored, etc) of the "matrix"; - template + template CudaLuDecomposition(const SparseMatrixPolicy& matrix) - : LuDecomposition(LuDecomposition::Create(matrix)) + : LuDecomposition(LuDecomposition::Create(matrix)) { /// Passing the class itself as an argument is not support by CUDA; /// Thus we generate a host struct first to save the pointers to @@ -75,22 +75,22 @@ namespace micm /// A is the sparse matrix to decompose /// L is the lower triangular matrix created by decomposition /// U is the upper triangular matrix created by decomposition - template typename SparseMatrixPolicy> - requires(CudaMatrix>&& VectorizableSparse>) void Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U) const; + template + requires(CudaMatrix&& VectorizableSparse) void Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U) const; - template typename SparseMatrixPolicy> - requires(!CudaMatrix>) void Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U) const; + template + requires(!CudaMatrix) void Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U) const; }; - template class SparseMatrixPolicy> - requires(CudaMatrix>&& VectorizableSparse>) void CudaLuDecomposition:: - Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const + template + requires(CudaMatrix&& VectorizableSparse) void CudaLuDecomposition:: + Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const { auto L_param = L.AsDeviceParam(); // we need to update lower matrix so it can't be constant and must be an lvalue auto U_param = U.AsDeviceParam(); // we need to update upper matrix so it can't be constant and must be an lvalue @@ -98,12 +98,12 @@ namespace micm } // call the function from the base class - template class SparseMatrixPolicy> - requires(!CudaMatrix>) void CudaLuDecomposition::Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U) const + template + requires(!CudaMatrix) void CudaLuDecomposition::Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U) const { - LuDecomposition::Decompose(A, L, U); + LuDecomposition::Decompose(A, L, U); } } // end of namespace micm diff --git a/include/micm/solver/cuda_rosenbrock.hpp b/include/micm/solver/cuda_rosenbrock.hpp index e72ce850b..23df355d4 100644 --- a/include/micm/solver/cuda_rosenbrock.hpp +++ b/include/micm/solver/cuda_rosenbrock.hpp @@ -35,9 +35,8 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, - class LinearSolverPolicy = CudaLinearSolver, + class LinearSolverPolicy = CudaLinearSolver, class ProcessSetPolicy = CudaProcessSet> class CudaRosenbrockSolver @@ -78,7 +77,7 @@ namespace micm const System& system, const std::vector processes, const RosenbrockSolverParameters& parameters, - const std::function, double)> create_linear_solver, + const std::function create_linear_solver, const std::function&, const std::map&)> create_process_set) : RosenbrockSolver( @@ -104,8 +103,8 @@ namespace micm micm::cuda::FreeConstData(this->devstruct_); }; - void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(CudaMatrix>&& VectorizableSparse>) + void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const + requires(CudaMatrix&& VectorizableSparse) { auto jacobian_param = jacobian.AsDeviceParam(); // we need to update jacobian so it can't be constant and must be an lvalue @@ -113,8 +112,8 @@ namespace micm } // call the function from the base class - void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(!CudaMatrix>) + void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const + requires(!CudaMatrix) { AlphaMinusJacobian(jacobian, alpha); } diff --git a/include/micm/solver/jit_linear_solver.hpp b/include/micm/solver/jit_linear_solver.hpp index cfd6a2274..0c5d55a56 100644 --- a/include/micm/solver/jit_linear_solver.hpp +++ b/include/micm/solver/jit_linear_solver.hpp @@ -19,10 +19,9 @@ namespace micm /// /// See LinearSolver class description for algorithm details /// The template parameter is the number of blocks (i.e. grid cells) in the block-diagonal matrix - template class SparseMatrixPolicy, class LuDecompositionPolicy = JitLuDecomposition> - class JitLinearSolver : public LinearSolver + template> + class JitLinearSolver : public LinearSolver { - std::shared_ptr compiler_; llvm::orc::ResourceTrackerSP solve_function_resource_tracker_; void (*solve_function_)(const double*, double*, const double*, const double*); @@ -39,8 +38,7 @@ namespace micm /// @param matrix Block-diagonal sparse matrix to create solver for /// @param initial_value Initial value for decomposed triangular matrix elements JitLinearSolver( - std::shared_ptr compiler, - const SparseMatrix>& matrix, + const SparseMatrixPolicy& matrix, double initial_value); ~JitLinearSolver(); @@ -49,25 +47,25 @@ namespace micm /// @param matrix Matrix that will be factored into lower and upper triangular matrices /// @param is_singular Flag that will be set to true if matrix is singular; false otherwise void Factor( - SparseMatrix>& matrix, - SparseMatrix>& lower_matrix, - SparseMatrix>& upper_matrix, + SparseMatrixPolicy& matrix, + SparseMatrixPolicy& lower_matrix, + SparseMatrixPolicy& upper_matrix, bool& is_singular); /// @brief Decompose the matrix into upper and lower triangular matrices and general JIT functions /// @param matrix Matrix that will be factored into lower and upper triangular matrices void Factor( - SparseMatrix>& matrix, - SparseMatrix>& lower_matrix, - SparseMatrix>& upper_matrix); + SparseMatrixPolicy& matrix, + SparseMatrixPolicy& lower_matrix, + SparseMatrixPolicy& upper_matrix); /// @brief Solve for x in Ax = b - template class MatrixPolicy> + template void Solve( - const MatrixPolicy& b, - MatrixPolicy& x, - SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix); + const MatrixPolicy& b, + MatrixPolicy& x, + SparseMatrixPolicy& lower_matrix, + SparseMatrixPolicy& upper_matrix); private: /// @brief Generates the JIT-ed Solve function diff --git a/include/micm/solver/jit_linear_solver.inl b/include/micm/solver/jit_linear_solver.inl index 3ce18247c..ff40277eb 100644 --- a/include/micm/solver/jit_linear_solver.inl +++ b/include/micm/solver/jit_linear_solver.inl @@ -5,39 +5,35 @@ namespace micm { - template class SparseMatrixPolicy, class LuDecompositionPolicy> + template inline JitLinearSolver::JitLinearSolver(JitLinearSolver &&other) - : LinearSolver(std::move(other)), - compiler_(std::move(other.compiler_)), + : LinearSolver(std::move(other)), solve_function_resource_tracker_(std::move(other.solve_function_resource_tracker_)), solve_function_(std::move(other.solve_function_)) { other.solve_function_ = NULL; } - template class SparseMatrixPolicy, class LuDecompositionPolicy> + template inline JitLinearSolver &JitLinearSolver::operator=(JitLinearSolver &&other) { - LinearSolver::operator=(std::move(other)); - compiler_ = std::move(other.compiler_); + LinearSolver::operator=(std::move(other)); solve_function_resource_tracker_ = std::move(other.solve_function_resource_tracker_); solve_function_ = std::move(other.solve_function_); other.solve_function_ = NULL; return *this; } - template class SparseMatrixPolicy, class LuDecompositionPolicy> + template inline JitLinearSolver::JitLinearSolver( - std::shared_ptr compiler, - const SparseMatrix> &matrix, + const SparseMatrixPolicy &matrix, double initial_value) - : LinearSolver( + : LinearSolver( matrix, initial_value, - [&](const SparseMatrixPolicy &m) -> LuDecompositionPolicy - { return LuDecompositionPolicy(compiler, m); }), - compiler_(compiler) + [&](const SparseMatrixPolicy &m) -> LuDecompositionPolicy + { return LuDecompositionPolicy(m); }) { solve_function_ = NULL; if (matrix.NumberOfBlocks() != L || matrix.GroupVectorSize() != L) @@ -51,7 +47,7 @@ namespace micm GenerateSolveFunction(); } - template class SparseMatrixPolicy, class LuDecompositionPolicy> + template inline JitLinearSolver::~JitLinearSolver() { if (solve_function_ != NULL) @@ -61,42 +57,42 @@ namespace micm } } - template class SparseMatrixPolicy, class LuDecompositionPolicy> + template inline void JitLinearSolver::Factor( - SparseMatrix> &matrix, - SparseMatrix> &lower_matrix, - SparseMatrix> &upper_matrix, + SparseMatrixPolicy &matrix, + SparseMatrixPolicy &lower_matrix, + SparseMatrixPolicy &upper_matrix, bool &is_singular) { - LinearSolver::Factor(matrix, lower_matrix, upper_matrix, is_singular); + LinearSolver::Factor(matrix, lower_matrix, upper_matrix, is_singular); } - template class SparseMatrixPolicy, class LuDecompositionPolicy> + template inline void JitLinearSolver::Factor( - SparseMatrix> &matrix, - SparseMatrix> &lower_matrix, - SparseMatrix> &upper_matrix) + SparseMatrixPolicy &matrix, + SparseMatrixPolicy &lower_matrix, + SparseMatrixPolicy &upper_matrix) { - LinearSolver::Factor(matrix, lower_matrix, upper_matrix); + LinearSolver::Factor(matrix, lower_matrix, upper_matrix); } - template class SparseMatrixPolicy, class LuDecompositionPolicy> - template class MatrixPolicy> + template + template inline void JitLinearSolver::Solve( - const MatrixPolicy &b, - MatrixPolicy &x, - SparseMatrixPolicy &lower_matrix, - SparseMatrixPolicy &upper_matrix) + const MatrixPolicy &b, + MatrixPolicy &x, + SparseMatrixPolicy &lower_matrix, + SparseMatrixPolicy &upper_matrix) { solve_function_( b.AsVector().data(), x.AsVector().data(), lower_matrix.AsVector().data(), upper_matrix.AsVector().data()); } - template class SparseMatrixPolicy, class LuDecompositionPolicy> + template inline void JitLinearSolver::GenerateSolveFunction() { std::string function_name = "linear_solve_" + GenerateRandomString(); - JitFunction func = JitFunction::Create(compiler_) + JitFunction func = JitFunction::Create() .SetName(function_name) .SetArguments({ { "b", JitType::DoublePtr }, { "x", JitType::DoublePtr }, @@ -104,10 +100,10 @@ namespace micm { "U", JitType::DoublePtr } }) .SetReturnType(JitType::Void); llvm::Type *double_type = func.GetType(JitType::Double); - auto Lij_yj = LinearSolver::Lij_yj_.begin(); - auto Uij_xj = LinearSolver::Uij_xj_.begin(); + auto Lij_yj = LinearSolver::Lij_yj_.begin(); + auto Uij_xj = LinearSolver::Uij_xj_.begin(); std::size_t offset = 0; - for (auto &nLij_Lii : LinearSolver::nLij_Lii_) + for (auto &nLij_Lii : LinearSolver::nLij_Lii_) { // the x vector is used for y values to conserve memory { @@ -157,8 +153,8 @@ namespace micm } offset += L; } - offset = L * LinearSolver::nUij_Uii_.size(); - for (auto &nUij_Uii : LinearSolver::nUij_Uii_) + offset = L * LinearSolver::nUij_Uii_.size(); + for (auto &nUij_Uii : LinearSolver::nUij_Uii_) { offset -= L; for (std::size_t i = 0; i < nUij_Uii.first; ++i) diff --git a/include/micm/solver/jit_lu_decomposition.hpp b/include/micm/solver/jit_lu_decomposition.hpp index 9bd40dc27..08a01e091 100644 --- a/include/micm/solver/jit_lu_decomposition.hpp +++ b/include/micm/solver/jit_lu_decomposition.hpp @@ -4,7 +4,6 @@ */ #pragma once -#include #include #include #include @@ -20,7 +19,6 @@ namespace micm template class JitLuDecomposition : public LuDecomposition { - std::shared_ptr compiler_; llvm::orc::ResourceTrackerSP decompose_function_resource_tracker_; void (*decompose_function_)(const double *, double *, double *); @@ -33,11 +31,8 @@ namespace micm JitLuDecomposition &operator=(JitLuDecomposition &&other); /// @brief Create a JITed LU decomposer for a given sparse matrix structure - /// @param compiler JIT compiler /// @param matrix Sparse matrix to create LU decomposer for - JitLuDecomposition( - std::shared_ptr compiler, - const SparseMatrix> &matrix); + JitLuDecomposition(const SparseMatrix> &matrix); ~JitLuDecomposition(); @@ -46,19 +41,19 @@ namespace micm /// @param lower The lower triangular matrix created by decomposition /// @param upper The upper triangular matrix created by decomposition /// @param is_singular Flag that will be set to true if A is singular; false otherwise - template class SparseMatrixPolicy> + template void Decompose( - const SparseMatrixPolicy &A, - SparseMatrixPolicy &lower, - SparseMatrixPolicy &upper, + const SparseMatrixPolicy &A, + SparseMatrixPolicy &lower, + SparseMatrixPolicy &upper, bool &is_singular) const; /// @brief Create sparse L and U matrices for a given A matrix /// @param A Sparse matrix that will be decomposed /// @param lower The lower triangular matrix created by decomposition /// @param upper The upper triangular matrix created by decomposition - template class SparseMatrixPolicy> - void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper) const; + template + void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper) const; private: /// @brief Generates a function to perform the LU decomposition for a specific matrix sparsity structure diff --git a/include/micm/solver/jit_lu_decomposition.inl b/include/micm/solver/jit_lu_decomposition.inl index 6b62a16d3..a1004c864 100644 --- a/include/micm/solver/jit_lu_decomposition.inl +++ b/include/micm/solver/jit_lu_decomposition.inl @@ -8,7 +8,6 @@ namespace micm template inline JitLuDecomposition::JitLuDecomposition(JitLuDecomposition &&other) : LuDecomposition(std::move(other)), - compiler_(std::move(other.compiler_)), decompose_function_resource_tracker_(std::move(other.decompose_function_resource_tracker_)), decompose_function_(std::move(other.decompose_function_)) { @@ -19,7 +18,6 @@ namespace micm inline JitLuDecomposition &JitLuDecomposition::operator=(JitLuDecomposition &&other) { LuDecomposition::operator=(std::move(other)); - compiler_ = std::move(other.compiler_); decompose_function_resource_tracker_ = std::move(other.decompose_function_resource_tracker_); decompose_function_ = std::move(other.decompose_function_); other.decompose_function_ = NULL; @@ -28,10 +26,8 @@ namespace micm template inline JitLuDecomposition::JitLuDecomposition( - std::shared_ptr compiler, const SparseMatrix> &matrix) - : LuDecomposition(LuDecomposition::Create>>(matrix)), - compiler_(compiler) + : LuDecomposition(LuDecomposition::Create>>(matrix)) { decompose_function_ = NULL; if (matrix.NumberOfBlocks() > L) @@ -59,7 +55,7 @@ namespace micm void JitLuDecomposition::GenerateDecomposeFunction() { std::string function_name = "lu_decompose_" + GenerateRandomString(); - JitFunction func = JitFunction::Create(compiler_) + JitFunction func = JitFunction::Create() .SetName(function_name) .SetArguments({ { "A matrix", JitType::DoublePtr }, { "lower matrix", JitType::DoublePtr }, @@ -190,22 +186,22 @@ namespace micm } template - template class SparseMatrixPolicy> + template void JitLuDecomposition::Decompose( - const SparseMatrixPolicy &A, - SparseMatrixPolicy &lower, - SparseMatrixPolicy &upper, + const SparseMatrixPolicy &A, + SparseMatrixPolicy &lower, + SparseMatrixPolicy &upper, bool &is_singular) const { - LuDecomposition::Decompose(A, lower, upper, is_singular); + LuDecomposition::Decompose(A, lower, upper, is_singular); } template - template class SparseMatrixPolicy> + template void JitLuDecomposition::Decompose( - const SparseMatrixPolicy &A, - SparseMatrixPolicy &lower, - SparseMatrixPolicy &upper) const + const SparseMatrixPolicy &A, + SparseMatrixPolicy &lower, + SparseMatrixPolicy &upper) const { decompose_function_(A.AsVector().data(), lower.AsVector().data(), upper.AsVector().data()); } diff --git a/include/micm/solver/jit_rosenbrock.hpp b/include/micm/solver/jit_rosenbrock.hpp index a9645c539..9f0a1ea0d 100644 --- a/include/micm/solver/jit_rosenbrock.hpp +++ b/include/micm/solver/jit_rosenbrock.hpp @@ -15,7 +15,6 @@ */ #pragma once -#include #include #include #include @@ -33,12 +32,11 @@ namespace micm template< template class MatrixPolicy = VectorMatrix, - template class SparseMatrixPolicy = VectorSparseMatrix, + class SparseMatrixPolicy = DefaultVectorSparseMatrix, class LinearSolverPolicy = JitLinearSolver, class ProcessSetPolicy = JitProcessSet> class JitRosenbrockSolver : public RosenbrockSolver { - std::shared_ptr compiler_; llvm::orc::ResourceTrackerSP function_resource_tracker_; using FuncPtr = void (*)(double*, const double); FuncPtr alpha_minus_jacobian_ = nullptr; @@ -48,7 +46,6 @@ namespace micm JitRosenbrockSolver& operator=(const JitRosenbrockSolver&) = delete; JitRosenbrockSolver(JitRosenbrockSolver&& other) : RosenbrockSolver(std::move(other)), - compiler_(std::move(other.compiler_)), function_resource_tracker_(std::move(other.function_resource_tracker_)), alpha_minus_jacobian_(std::move(other.alpha_minus_jacobian_)) { @@ -58,7 +55,6 @@ namespace micm JitRosenbrockSolver& operator=(JitRosenbrockSolver&& other) { RosenbrockSolver::operator=(std::move(other)); - compiler_ = std::move(other.compiler_); function_resource_tracker_ = std::move(other.function_resource_tracker_); alpha_minus_jacobian_ = std::move(other.alpha_minus_jacobian_); other.alpha_minus_jacobian_ = NULL; @@ -69,7 +65,6 @@ namespace micm /// @param system The chemical system to create the solver for /// @param processes The collection of chemical processes that will be applied during solving JitRosenbrockSolver( - std::shared_ptr compiler, const System& system, const std::vector& processes, const RosenbrockSolverParameters& parameters) @@ -77,14 +72,13 @@ namespace micm system, processes, parameters, - [&](const SparseMatrixPolicy& matrix, double initial_value) -> LinearSolverPolicy { - return LinearSolverPolicy{ compiler, matrix, initial_value }; + [&](const SparseMatrixPolicy& matrix, double initial_value) -> LinearSolverPolicy { + return LinearSolverPolicy{ matrix, initial_value }; }, [&](const std::vector& processes, const std::map& variable_map) -> ProcessSetPolicy { - return ProcessSetPolicy{ compiler, processes, variable_map }; - }), - compiler_(compiler) + return ProcessSetPolicy{ processes, variable_map }; + }) { MatrixPolicy temp{}; if (temp.GroupVectorSize() != parameters.number_of_grid_cells_) @@ -110,7 +104,7 @@ namespace micm /// @brief compute [alpha * I - dforce_dy] /// @param jacobian Jacobian matrix (dforce_dy) /// @param alpha - void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const + void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const { double a = alpha; if (alpha_minus_jacobian_) @@ -133,7 +127,7 @@ namespace micm std::size_t number_of_nonzero_jacobian_elements = jacobian.AsVector().size(); // Create the JitFunction with the modified name - JitFunction func = JitFunction::Create(compiler_) + JitFunction func = JitFunction::Create() .SetName("alpha_minus_jacobian") .SetArguments({ { "jacobian", JitType::DoublePtr }, { "alpha", JitType::Double } }) .SetReturnType(JitType::Void); diff --git a/include/micm/solver/linear_solver.hpp b/include/micm/solver/linear_solver.hpp index 60d9cabbd..6bf5e9aaa 100644 --- a/include/micm/solver/linear_solver.hpp +++ b/include/micm/solver/linear_solver.hpp @@ -24,7 +24,7 @@ namespace micm /// @brief A general-use block-diagonal sparse-matrix linear solver /// /// The sparsity pattern of each block in the block diagonal matrix is the same. - template class SparseMatrixPolicy, class LuDecompositionPolicy = LuDecomposition> + template class LinearSolver { protected: @@ -60,43 +60,42 @@ namespace micm /// @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 - LinearSolver(const SparseMatrixPolicy& matrix, T initial_value); + LinearSolver(const SparseMatrixPolicy& matrix, typename SparseMatrixPolicy::value_type initial_value); /// @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 /// @param create_lu_decomp Function to create an LU Decomposition object that adheres to LuDecompositionPolicy LinearSolver( - const SparseMatrixPolicy& matrix, - T initial_value, - const std::function&)> create_lu_decomp); + const SparseMatrixPolicy& matrix, + typename SparseMatrixPolicy::value_type initial_value, + const std::function create_lu_decomp); /// @brief Decompose the matrix into upper and lower triangular matrices - void - Factor(const SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix); + void Factor(const SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix); /// @brief Decompose the matrix into upper and lower triangular matrices /// @param matrix Matrix to decompose into lower and upper triangular matrices /// @param is_singular Flag that is set to true if matrix is singular; false otherwise void Factor( - const SparseMatrixPolicy& matrix, - SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix, + const SparseMatrixPolicy& matrix, + SparseMatrixPolicy& lower_matrix, + SparseMatrixPolicy& upper_matrix, bool& is_singular); /// @brief Solve for x in Ax = b - template class MatrixPolicy> - requires(!VectorizableDense> || !VectorizableSparse>) void Solve( - const MatrixPolicy& b, - MatrixPolicy& x, - const SparseMatrixPolicy& lower_matrix, - const SparseMatrixPolicy& upper_matrix) const; - template class MatrixPolicy> - requires(VectorizableDense>&& VectorizableSparse>) void Solve( - const MatrixPolicy& b, - MatrixPolicy& x, - const SparseMatrixPolicy& lower_matrix, - const SparseMatrixPolicy& upper_matrix) const; + template + requires(!VectorizableDense || !VectorizableSparse) void Solve( + const MatrixPolicy& b, + MatrixPolicy& x, + const SparseMatrixPolicy& lower_matrix, + const SparseMatrixPolicy& upper_matrix) const; + template + requires(VectorizableDense&& VectorizableSparse) void Solve( + const MatrixPolicy& b, + MatrixPolicy& x, + const SparseMatrixPolicy& lower_matrix, + const SparseMatrixPolicy& upper_matrix) const; }; } // namespace micm diff --git a/include/micm/solver/linear_solver.inl b/include/micm/solver/linear_solver.inl index 2128252f9..d252b9aa6 100644 --- a/include/micm/solver/linear_solver.inl +++ b/include/micm/solver/linear_solver.inl @@ -2,17 +2,18 @@ * * SPDX-License-Identifier: Apache-2.0 */ + namespace micm { - template class MatrixPolicy> - inline std::vector DiagonalMarkowitzReorder(const MatrixPolicy& matrix) + template + inline std::vector DiagonalMarkowitzReorder(const MatrixPolicy& matrix) { MICM_PROFILE_FUNCTION(); const std::size_t order = matrix.NumRows(); std::vector perm(order); for (std::size_t i = 0; i < order; ++i) perm[i] = i; - MatrixPolicy pattern = matrix; + MatrixPolicy pattern = matrix; for (std::size_t row = 0; row < (order - 1); ++row) { std::size_t beta = std::pow((order - 1), 2); @@ -50,23 +51,23 @@ namespace micm return perm; } - template class SparseMatrixPolicy, class LuDecompositionPolicy> - inline LinearSolver::LinearSolver( - const SparseMatrixPolicy& matrix, - T initial_value) - : LinearSolver( + template + inline LinearSolver::LinearSolver( + const SparseMatrixPolicy& matrix, + typename SparseMatrixPolicy::value_type initial_value) + : LinearSolver( matrix, initial_value, - [](const SparseMatrixPolicy& m) -> LuDecompositionPolicy - { return LuDecomposition::Create(m); }) + [](const SparseMatrixPolicy& m) -> LuDecompositionPolicy + { return LuDecomposition::Create(m); }) { } - template class SparseMatrixPolicy, class LuDecompositionPolicy> - inline LinearSolver::LinearSolver( - const SparseMatrixPolicy& matrix, - T initial_value, - const std::function&)> create_lu_decomp) + template + inline LinearSolver::LinearSolver( + const SparseMatrixPolicy& matrix, + typename SparseMatrixPolicy::value_type initial_value, + const std::function create_lu_decomp) : nLij_Lii_(), Lij_yj_(), nUij_Uii_(), @@ -75,7 +76,7 @@ namespace micm { MICM_PROFILE_FUNCTION(); - auto lu = lu_decomp_.template GetLUMatrices(matrix, initial_value); + auto lu = lu_decomp_.template GetLUMatrices(matrix, initial_value); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); for (std::size_t i = 0; i < lower_matrix.NumRows(); ++i) @@ -108,39 +109,38 @@ namespace micm } }; - template class SparseMatrixPolicy, class LuDecompositionPolicy> - inline void LinearSolver::Factor( - const SparseMatrixPolicy& matrix, - SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix) + template + inline void LinearSolver::Factor( + const SparseMatrixPolicy& matrix, + SparseMatrixPolicy& lower_matrix, + SparseMatrixPolicy& upper_matrix) { MICM_PROFILE_FUNCTION(); - lu_decomp_.template Decompose(matrix, lower_matrix, upper_matrix); + lu_decomp_.template Decompose(matrix, lower_matrix, upper_matrix); } - template class SparseMatrixPolicy, class LuDecompositionPolicy> - inline void LinearSolver::Factor( - const SparseMatrixPolicy& matrix, - SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix, + template + inline void LinearSolver::Factor( + const SparseMatrixPolicy& matrix, + SparseMatrixPolicy& lower_matrix, + SparseMatrixPolicy& upper_matrix, bool& is_singular) { MICM_PROFILE_FUNCTION(); - lu_decomp_.template Decompose(matrix, lower_matrix, upper_matrix, is_singular); + lu_decomp_.template Decompose(matrix, lower_matrix, upper_matrix, is_singular); } - template class SparseMatrixPolicy, class LuDecompositionPolicy> - template class MatrixPolicy> + template + template requires( - !VectorizableDense> || - !VectorizableSparse>) inline void LinearSolver:: - Solve( - const MatrixPolicy& b, - MatrixPolicy& x, - const SparseMatrixPolicy& lower_matrix, - const SparseMatrixPolicy& upper_matrix) const + !VectorizableDense || + !VectorizableSparse) inline void LinearSolver::Solve( + const MatrixPolicy& b, + MatrixPolicy& x, + const SparseMatrixPolicy& lower_matrix, + const SparseMatrixPolicy& upper_matrix) const { MICM_PROFILE_FUNCTION(); @@ -189,15 +189,15 @@ namespace micm } } - template class SparseMatrixPolicy, class LuDecompositionPolicy> - template class MatrixPolicy> - requires(VectorizableDense>&& VectorizableSparse< - SparseMatrixPolicy>) inline void LinearSolver:: + template + template + requires(VectorizableDense&& VectorizableSparse< + SparseMatrixPolicy>) inline void LinearSolver:: Solve( - const MatrixPolicy& b, - MatrixPolicy& x, - const SparseMatrixPolicy& lower_matrix, - const SparseMatrixPolicy& upper_matrix) const + const MatrixPolicy& b, + MatrixPolicy& x, + const SparseMatrixPolicy& lower_matrix, + const SparseMatrixPolicy& upper_matrix) const { MICM_PROFILE_FUNCTION(); const std::size_t n_cells = b.GroupVectorSize(); diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp index 9e8906ba2..3173f1053 100644 --- a/include/micm/solver/lu_decomposition.hpp +++ b/include/micm/solver/lu_decomposition.hpp @@ -78,56 +78,50 @@ namespace micm /// @brief Construct an LU decomposition algorithm for a given sparse matrix /// @param matrix Sparse matrix - template - LuDecomposition(const SparseMatrix& matrix); + template + LuDecomposition(const SparseMatrixPolicy& matrix); /// @brief Create an LU decomposition algorithm for a given sparse matrix policy /// @param matrix Sparse matrix - template class SparseMatrixPolicy> - static LuDecomposition Create(const SparseMatrixPolicy& matrix); - template + template 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 class SparseMatrixPolicy> - static std::pair, SparseMatrixPolicy> GetLUMatrices( - const SparseMatrixPolicy& A, - T initial_value); - template - static std::pair GetLUMatrices(const SparseMatrixPolicy& A, T initial_value); + template + static std::pair GetLUMatrices(const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value); /// @brief Perform an LU decomposition on a given A matrix /// @param A Sparse matrix to decompose /// @param L The lower triangular matrix created by decomposition /// @param U The upper triangular matrix created by decomposition - template class SparseMatrixPolicy> - void Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const; + template + void Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const; /// @brief Perform an LU decomposition on a given A matrix /// @param A Sparse matrix to decompose /// @param L The lower triangular matrix created by decomposition /// @param U The upper triangular matrix created by decomposition /// @param is_singular Flag that is set to true if A is singular; false otherwise - template class SparseMatrixPolicy> - requires(!VectorizableSparse>) void Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U, + template + requires(!VectorizableSparse) void Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U, bool& is_singular) const; - template class SparseMatrixPolicy> - requires(VectorizableSparse>) void Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U, + template + requires(VectorizableSparse) void Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U, bool& is_singular) const; private: /// @brief Initialize arrays for the LU decomposition /// @param A Sparse matrix to decompose - template - void Initialize(const SparseMatrixPolicy& matrix, T initial_value); + template + void Initialize(const SparseMatrixPolicy& matrix, auto initial_value); }; } // namespace micm diff --git a/include/micm/solver/lu_decomposition.inl b/include/micm/solver/lu_decomposition.inl index 70cf2f915..d9803d42b 100644 --- a/include/micm/solver/lu_decomposition.inl +++ b/include/micm/solver/lu_decomposition.inl @@ -9,35 +9,27 @@ namespace micm { } - template - inline LuDecomposition::LuDecomposition(const SparseMatrix& matrix) + template + inline LuDecomposition::LuDecomposition(const SparseMatrixPolicy& matrix) { - Initialize(matrix); + Initialize(matrix, typename SparseMatrixPolicy::value_type()); } - template class SparseMatrixPolicy> - inline LuDecomposition LuDecomposition::Create(const SparseMatrixPolicy& matrix) - { - LuDecomposition lu_decomp{}; - lu_decomp.Initialize>(matrix, T{}); - return lu_decomp; - } - - template + template inline LuDecomposition LuDecomposition::Create(const SparseMatrixPolicy& matrix) { LuDecomposition lu_decomp{}; - lu_decomp.Initialize(matrix, T{}); + lu_decomp.Initialize(matrix, typename SparseMatrixPolicy::value_type()); return lu_decomp; } - template - inline void LuDecomposition::Initialize(const SparseMatrixPolicy& matrix, T initial_value) + template + inline void LuDecomposition::Initialize(const SparseMatrixPolicy& matrix, auto initial_value) { MICM_PROFILE_FUNCTION(); std::size_t n = matrix.NumRows(); - auto LU = GetLUMatrices(matrix, initial_value); + auto LU = GetLUMatrices(matrix, initial_value); const auto& L_row_start = LU.first.RowStartVector(); const auto& L_row_ids = LU.first.RowIdsVector(); const auto& U_row_start = LU.second.RowStartVector(); @@ -107,18 +99,10 @@ namespace micm } } - template class SparseMatrixPolicy> - inline std::pair, SparseMatrixPolicy> LuDecomposition::GetLUMatrices( - const SparseMatrixPolicy& A, - T initial_value) - { - return GetLUMatrices>(A, initial_value); - } - - template + template inline std::pair LuDecomposition::GetLUMatrices( const SparseMatrixPolicy& A, - T initial_value) + typename SparseMatrixPolicy::value_type initial_value) { MICM_PROFILE_FUNCTION(); @@ -177,19 +161,19 @@ namespace micm return LU; } - template class SparseMatrixPolicy> - inline void LuDecomposition::Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) + template + inline void LuDecomposition::Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const { bool is_singular = false; - Decompose(A, L, U, is_singular); + Decompose(A, L, U, is_singular); } - template class SparseMatrixPolicy> - requires(!VectorizableSparse>) inline void LuDecomposition::Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U, + template + requires(!VectorizableSparse) inline void LuDecomposition::Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U, bool& is_singular) const { MICM_PROFILE_FUNCTION(); @@ -248,11 +232,11 @@ namespace micm } } - template class SparseMatrixPolicy> - requires(VectorizableSparse>) inline void LuDecomposition::Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U, + template + requires(VectorizableSparse) inline void LuDecomposition::Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U, bool& is_singular) const { MICM_PROFILE_FUNCTION(); diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index cdfe2bba3..4a73b2e5b 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -92,8 +92,8 @@ namespace micm /// The template parameter is the type of matrix to use template< template class MatrixPolicy = Matrix, - template class SparseMatrixPolicy = StandardSparseMatrix, - class LinearSolverPolicy = LinearSolver, + class SparseMatrixPolicy = StandardSparseMatrix, + class LinearSolverPolicy = LinearSolver, class ProcessSetPolicy = ProcessSet> class RosenbrockSolver { @@ -141,7 +141,7 @@ namespace micm const System& system, const std::vector& processes, const RosenbrockSolverParameters& parameters, - const std::function, double)> create_linear_solver, + const std::function create_linear_solver, const std::function&, const std::map&)> create_process_set); @@ -149,12 +149,12 @@ namespace micm /// @brief Returns a state object for use with the solver /// @return A object that can hold the full state of the chemical system - virtual State GetState() const; + virtual State, SparseMatrixPolicy> GetState() const; /// @brief Advances the given step over the specified time step /// @param time_step Time [s] to advance the state by /// @return A struct containing results and a status code - SolverResult Solve(double time_step, State& state) noexcept; + SolverResult Solve(double time_step, State, SparseMatrixPolicy>& state) noexcept; /// @brief Calculate a chemical forcing /// @param rate_constants List of rate constants for each needed species @@ -168,14 +168,14 @@ namespace micm /// @brief compute [alpha * I - dforce_dy] /// @param jacobian Jacobian matrix (dforce_dy) /// @param alpha - void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(!VectorizableSparse>); - void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(VectorizableSparse>); + void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const + requires(!VectorizableSparse); + void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const + requires(VectorizableSparse); /// @brief Update the rate constants for the environment state /// @param state The current state of the chemical system - void UpdateState(State& state); + void UpdateState(State, SparseMatrixPolicy>& state); /// @brief Compute the derivative of the forcing w.r.t. each chemical, and return the negative jacobian /// @param rate_constants List of rate constants for each needed species @@ -184,7 +184,7 @@ namespace micm virtual void CalculateNegativeJacobian( const MatrixPolicy& rate_constants, const MatrixPolicy& number_densities, - SparseMatrixPolicy& jacobian); + SparseMatrixPolicy& jacobian); /// @brief Prepare the linear solver /// @param H time step (seconds) @@ -198,7 +198,7 @@ namespace micm bool& singular, const MatrixPolicy& number_densities, SolverStats& stats, - State& state); + State, SparseMatrixPolicy>& state); /// @brief Computes the scaled norm of the vector errors /// @param y the original vector diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 5f84f4870..9a48b3c07 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -78,7 +78,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -94,7 +93,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -106,7 +104,7 @@ namespace micm system, processes, parameters, - [](const SparseMatrixPolicy& matrix, double initial_value) -> LinearSolverPolicy { + [](const SparseMatrixPolicy& matrix, double initial_value) -> LinearSolverPolicy { return LinearSolverPolicy{ matrix, initial_value }; }, [](const std::vector& processes, @@ -119,7 +117,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -127,7 +124,7 @@ namespace micm const System& system, const std::vector& processes, const RosenbrockSolverParameters& parameters, - const std::function, double)> create_linear_solver, + const std::function create_linear_solver, const std::function&, const std::map&)> create_process_set) : processes_(processes), @@ -172,7 +169,7 @@ namespace micm MatrixPolicy unsorted_jac_non_zeros(system.StateSize(), system.StateSize(), 0); for (auto& elem : unsorted_jac_elements) unsorted_jac_non_zeros[elem.first][elem.second] = 1; - auto reorder_map = DiagonalMarkowitzReorder(unsorted_jac_non_zeros); + auto reorder_map = DiagonalMarkowitzReorder>(unsorted_jac_non_zeros); state_reordering = [=](const std::vector& variables, const std::size_t i) { return variables[reorder_map[i]]; }; @@ -237,27 +234,25 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> - inline State + inline State, SparseMatrixPolicy> RosenbrockSolver::GetState() const { - return State{ state_parameters_ }; + return State, SparseMatrixPolicy>{ state_parameters_ }; } template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline typename RosenbrockSolver::SolverResult RosenbrockSolver::Solve( double time_step, - State& state) noexcept + State, SparseMatrixPolicy>& state) noexcept { MICM_PROFILE_FUNCTION(); @@ -367,7 +362,7 @@ namespace micm K[stage].Axpy(parameters_.c_[stage_combinations + j] / H, K[j]); } temp.Copy(K[stage]); - linear_solver_.template Solve(temp, K[stage], state.lower_matrix_, state.upper_matrix_); + linear_solver_.template Solve>(temp, K[stage], state.lower_matrix_, state.upper_matrix_); stats.solves_ += 1; } @@ -455,7 +450,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -466,19 +460,18 @@ namespace micm { MICM_PROFILE_FUNCTION(); std::fill(forcing.AsVector().begin(), forcing.AsVector().end(), 0.0); - process_set_.template AddForcingTerms(rate_constants, number_densities, forcing); + process_set_.template AddForcingTerms>(rate_constants, number_densities, forcing); } template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::AlphaMinusJacobian( - SparseMatrixPolicy& jacobian, - const double& alpha) const requires(!VectorizableSparse>) + SparseMatrixPolicy& jacobian, + const double& alpha) const requires(!VectorizableSparse) { MICM_PROFILE_FUNCTION(); @@ -493,13 +486,12 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::AlphaMinusJacobian( - SparseMatrixPolicy& jacobian, - const double& alpha) const requires(VectorizableSparse>) + SparseMatrixPolicy& jacobian, + const double& alpha) const requires(VectorizableSparse) { MICM_PROFILE_FUNCTION(); @@ -516,7 +508,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -524,7 +515,7 @@ namespace micm RosenbrockSolver::CalculateNegativeJacobian( const MatrixPolicy& rate_constants, const MatrixPolicy& number_densities, - SparseMatrixPolicy& jacobian) + SparseMatrixPolicy& jacobian) { MICM_PROFILE_FUNCTION(); @@ -535,12 +526,11 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::UpdateState( - State& state) + State, SparseMatrixPolicy>& state) { Process::UpdateState(processes_, state); } @@ -548,7 +538,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -558,7 +547,7 @@ namespace micm bool& singular, const MatrixPolicy& number_densities, SolverStats& stats, - State& state) + State, SparseMatrixPolicy>& state) { MICM_PROFILE_FUNCTION(); @@ -592,7 +581,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -632,7 +620,6 @@ namespace micm template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 941eadf86..46fc79c12 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -17,9 +17,9 @@ namespace micm /// @brief Rosenbrock solver parameters struct RosenbrockSolverParameters { - size_t stages_{}; - size_t upper_limit_tolerance_{}; - size_t max_number_of_steps_{ 1000 }; + std::size_t stages_{}; + std::size_t upper_limit_tolerance_{}; + std::size_t max_number_of_steps_{ 1000 }; double round_off_{ std::numeric_limits::epsilon() }; // Unit roundoff (1+round_off)>1 double factor_min_{ 0.2 }; // solver step size minimum boundary @@ -59,7 +59,7 @@ namespace micm std::vector absolute_tolerance_; double relative_tolerance_{ 1e-4 }; - size_t number_of_grid_cells_{ 1 }; // Number of grid cells to solve simultaneously + std::size_t number_of_grid_cells_{ 1 }; // Number of grid cells to solve simultaneously bool reorder_state_{ true }; // Reorder state during solver construction to minimize LU fill-in bool check_singularity_{ false }; // Check for singular A matrix in linear solve of A x = b bool ignore_unused_species_{ false }; // Allow unused species to be included in state and solve @@ -72,35 +72,35 @@ namespace micm /// @param reorder_state /// @return static RosenbrockSolverParameters TwoStageRosenbrockParameters( - size_t number_of_grid_cells = 1, + std::size_t number_of_grid_cells = 1, bool reorder_state = true); /// @brief an L-stable method, 3 stages, order 3, 2 function evaluations /// @param number_of_grid_cells /// @param reorder_state /// @return static RosenbrockSolverParameters ThreeStageRosenbrockParameters( - size_t number_of_grid_cells = 1, + std::size_t number_of_grid_cells = 1, bool reorder_state = true); /// @brief L-stable rosenbrock method of order 4, with 4 stages /// @param number_of_grid_cells /// @param reorder_state /// @return static RosenbrockSolverParameters FourStageRosenbrockParameters( - size_t number_of_grid_cells = 1, + std::size_t number_of_grid_cells = 1, bool reorder_state = true); /// @brief A stiffly-stable method, 4 stages, order 3 /// @param number_of_grid_cells /// @param reorder_state /// @return static RosenbrockSolverParameters FourStageDifferentialAlgebraicRosenbrockParameters( - size_t number_of_grid_cells = 1, + std::size_t number_of_grid_cells = 1, bool reorder_state = true); /// @brief stiffly-stable rosenbrock method of order 4, with 6 stages /// @param number_of_grid_cells /// @param reorder_state /// @return static RosenbrockSolverParameters SixStageDifferentialAlgebraicRosenbrockParameters( - size_t number_of_grid_cells = 1, + std::size_t number_of_grid_cells = 1, bool reorder_state = true); private: @@ -160,7 +160,7 @@ namespace micm } inline RosenbrockSolverParameters RosenbrockSolverParameters::TwoStageRosenbrockParameters( - size_t number_of_grid_cells, + std::size_t number_of_grid_cells, bool reorder_state) { // an L-stable method, 2 stages, order 2 @@ -204,7 +204,7 @@ namespace micm } inline RosenbrockSolverParameters RosenbrockSolverParameters::ThreeStageRosenbrockParameters( - size_t number_of_grid_cells, + std::size_t number_of_grid_cells, bool reorder_state) { // an L-stable method, 3 stages, order 3, 2 function evaluations @@ -260,7 +260,7 @@ namespace micm } inline RosenbrockSolverParameters RosenbrockSolverParameters::FourStageRosenbrockParameters( - size_t number_of_grid_cells, + std::size_t number_of_grid_cells, bool reorder_state) { // L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES @@ -329,7 +329,7 @@ namespace micm } inline RosenbrockSolverParameters RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters( - size_t number_of_grid_cells, + std::size_t number_of_grid_cells, bool reorder_state) { // A STIFFLY-STABLE METHOD, 4 stages, order 3 @@ -384,7 +384,7 @@ namespace micm } inline RosenbrockSolverParameters RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters( - size_t number_of_grid_cells, + std::size_t number_of_grid_cells, bool reorder_state) { // STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES diff --git a/include/micm/solver/solver.hpp b/include/micm/solver/solver.hpp new file mode 100644 index 000000000..33dea939c --- /dev/null +++ b/include/micm/solver/solver.hpp @@ -0,0 +1,68 @@ +/* Copyright (C) 2023-2024 National Center for Atmospheric Research + * + * SPDX-License-Identifier: Apache-2.0 + */ +#pragma once + + +namespace micm +{ + + template + class Solver + { + private: + std::size_t number_of_grid_cells_; + std::size_t number_of_species_; + std::size_t number_of_reactions_; + StateParameters state_parameters_; + SolverPolicy solver_; + + public: + + 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), + number_of_grid_cells_(number_of_grid_cells), + number_of_species_(number_of_species), + number_of_reactions_(number_of_reactions), + state_parameters_(state_parameters) + { + } + + void Solve(double time_step, StatePolicy& state) + { + solver_.Solve(time_step, state); + } + + /// @brief Returns the number of grid cells + /// @return + std::size_t GetNumberOfGridCells() const + { + return number_of_grid_cells_; + } + + /// @brief Returns the number of species + /// @return + std::size_t GetNumberOfSpecies() const + { + return number_of_species_; + } + + /// @brief Returns the number of reactions + /// @return + std::size_t GetNumberOfReactions() const + { + return number_of_reactions_; + } + + StatePolicy GetState() const + { + return StatePolicy(state_parameters_); + } + }; +} // namespace micm diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp new file mode 100644 index 000000000..5b91abb69 --- /dev/null +++ b/include/micm/solver/solver_builder.hpp @@ -0,0 +1,116 @@ +/* Copyright (C) 2023-2024 National Center for Atmospheric Research + * + * SPDX-License-Identifier: Apache-2.0 + */ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace micm +{ + using SolverParameters = std::variant; + + template + class SolverBuilder + { + protected: + System system_; + std::size_t number_of_grid_cells_; + std::vector reactions_; + SolverParameters options_; + bool ignore_unused_species_ = true; + bool reorder_state_ = true; + + public: + /// @brief Set the chemical system + /// @param system + /// @return + SolverBuilder& SetSystem(const System& system); + + /// @brief Set the reactions + /// @param reactions + /// @return + SolverBuilder& SetReactions(const std::vector& reactions); + + /// @brief Set the number of grid cells + /// @param number_of_grid_cells + /// @return + SolverBuilder& SetNumberOfGridCells(int number_of_grid_cells); + + /// @brief Choose a rosenbrock solver + /// @param options + /// @return + SolverBuilder& SetSolverParameters(const RosenbrockSolverParameters& options); + + /// @brief Choose a backward euler solver + /// @param options + /// @return + SolverBuilder& SetSolverParameters(const BackwardEulerSolverParameters& options); + + /// @brief + /// @return + auto Build(); + + protected: + /// @brief + /// @return + virtual Solver, ProcessSet>, State> BuildBackwardEulerSolver() = 0; + + virtual ~SolverBuilder() = default; + + /// @brief + /// @tparam ProcessSetPolicy + /// @return + template + void UnusedSpeciesCheck(); + + /// @brief Get a species map properly ordered + /// @return + template + std::map GetSpeciesMap() const; + + /// @brief Set the absolute tolerances per species + /// @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 + std::vector GetCustomParameterLabels() const; + + std::vector GetJacobianDiagonalElements(auto jacobian) const; + }; + + template + class CpuSolverBuilder : public SolverBuilder + { + public: + Solver, ProcessSet>, State> BuildBackwardEulerSolver() override; + }; + + template + class GpuSolverBuilder : public SolverBuilder + { + public: + Solver, ProcessSet>, State> BuildBackwardEulerSolver() override; + }; + +} // namespace micm + +#include "solver_builder.inl" \ No newline at end of file diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl new file mode 100644 index 000000000..504d8d9d1 --- /dev/null +++ b/include/micm/solver/solver_builder.inl @@ -0,0 +1,277 @@ +/* Copyright (C) 2023-2024 National Center for Atmospheric Research + * + * SPDX-License-Identifier: Apache-2.0 + */ + +enum class MicmSolverBuilderErrc +{ + UnusedSpecies = 1, // Unused species present in the chemical system +}; + +namespace std +{ + template<> + struct is_error_condition_enum : true_type + { + }; +} // namespace std + +namespace +{ + class SolverBuilderErrorCategory : public std::error_category + { + public: + const char* name() const noexcept override + { + return "MICM Solver Builder"; + } + std::string message(int ev) const override + { + switch (static_cast(ev)) + { + case MicmSolverBuilderErrc::UnusedSpecies: + return "Unused species present in the chemical system. Use the ignore_unused_species_ parameter to allow unused " + "species in the solve."; + default: return "Unknown error"; + } + } + }; + + const SolverBuilderErrorCategory solverBuilderErrorCategory{}; +} // namespace + +inline std::error_code make_error_code(MicmSolverBuilderErrc e) +{ + return { static_cast(e), solverBuilderErrorCategory }; +} + +namespace micm +{ + template + inline SolverBuilder& + SolverBuilder::SetSystem(const System& system) + { + system_ = system; + return *this; + } + + template + inline SolverBuilder& + SolverBuilder::SetReactions(const std::vector& reactions) + { + reactions_ = reactions; + return *this; + } + + template + inline SolverBuilder& + SolverBuilder::SetNumberOfGridCells(int number_of_grid_cells) + { + number_of_grid_cells_ = number_of_grid_cells; + return *this; + } + + template + inline SolverBuilder& + SolverBuilder::SetSolverParameters(const RosenbrockSolverParameters& options) + { + if (!std::holds_alternative(options_)) + { + throw std::runtime_error("Solver type already set"); + } + options_.emplace(options); + return *this; + } + + template + inline SolverBuilder& + SolverBuilder::SetSolverParameters(const BackwardEulerSolverParameters& options) + { + if (!std::holds_alternative(options_)) + { + throw std::runtime_error("Solver type already set"); + } + options_.emplace(options); + return *this; + } + + template + inline auto SolverBuilder::Build() + { + if (std::holds_alternative(options_)) + { + throw std::runtime_error("Not implemented yet"); + } + if (std::holds_alternative(options_)) + { + return BuildBackwardEulerSolver(); + } + + throw std::runtime_error("No solver type set"); + } + + template + template + inline void SolverBuilder::UnusedSpeciesCheck() + { + if (ignore_unused_species_) + { + return; + } + + auto used_species = ProcessSetPolicy::SpeciesUsed(reactions_); + auto available_species = system_.UniqueNames(); + std::sort(available_species.begin(), available_species.end()); + std::set unused_species; + std::set_difference( + available_species.begin(), + available_species.end(), + used_species.begin(), + used_species.end(), + std::inserter(unused_species, unused_species.begin())); + if (unused_species.size() > 0) + { + std::string err_msg = "Unused species in chemical system:"; + for (auto& species : unused_species) + err_msg += " '" + species + "'"; + err_msg += "."; + throw std::system_error(make_error_code(MicmSolverBuilderErrc::UnusedSpecies), err_msg); + } + } + + template + template + inline std::map SolverBuilder::GetSpeciesMap() const + { + std::map species_map; + std::function& variables, const std::size_t i)> state_reordering; + std::size_t index = 0; + for (auto& name : system_.UniqueNames()) + species_map[name] = index++; + + if (reorder_state_) + { + // get unsorted Jacobian non-zero elements + auto unsorted_process_set = ProcessSetPolicy(reactions_, species_map); + auto unsorted_jac_elements = unsorted_process_set.NonZeroJacobianElements(); + + using Matrix = typename DenseMatrixPolicy::IntMatrix; + Matrix unsorted_jac_non_zeros(system_.StateSize(), system_.StateSize(), 0); + for (auto& elem : unsorted_jac_elements) + unsorted_jac_non_zeros[elem.first][elem.second] = 1; + auto reorder_map = DiagonalMarkowitzReorder(unsorted_jac_non_zeros); + + state_reordering = [=](const std::vector& variables, const std::size_t i) + { return variables[reorder_map[i]]; }; + + index = 0; + for (auto& name : system_.UniqueNames(state_reordering)) + species_map[name] = index++; + } + + return species_map; + } + + template + inline void SolverBuilder::SetAbsoluteTolerances( + std::vector& tolerances, + const std::map& species_map) const + { + // if the tolerances aren't already set, initialize them and then set based off of information in the system + if (tolerances.size() != species_map.size()) + { + tolerances = std::vector(species_map.size(), 1e-3); + for (auto& species : system_.gas_phase_.species_) + { + if (species.HasProperty("absolute tolerance")) + { + tolerances[species_map.at(species.name_)] = species.template GetProperty("absolute tolerance"); + } + } + for (auto& phase : system_.phases_) + { + for (auto& species : phase.second.species_) + { + if (species.HasProperty("absolute tolerance")) + { + tolerances[species_map.at(species.name_)] = species.template GetProperty("absolute tolerance"); + } + } + } + } + } + + template + inline std::vector SolverBuilder::GetCustomParameterLabels() const + { + std::vector 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; + } + + 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 Solver< + BackwardEuler, ProcessSet>, + State> + CpuSolverBuilder::BuildBackwardEulerSolver() + { + using ProcessSetPolicy = ProcessSet; + using LinearSolverPolicy = LinearSolver; + + auto parameters = std::get(this->options_); + auto species_map = this->template GetSpeciesMap(); + auto labels = this->GetCustomParameterLabels(); + std::size_t number_of_species = this->system_.StateSize(); + + this->template UnusedSpeciesCheck(); + this->SetAbsoluteTolerances(parameters.absolute_tolerance_, species_map); + + ProcessSetPolicy process_set(this->reactions_, species_map); + auto diagonal_elements = process_set.NonZeroJacobianElements(); + auto jacobian = BuildJacobian(diagonal_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); + + std::vector variable_names{ number_of_species }; + for (auto& species_pair : species_map) + variable_names[species_pair.second] = species_pair.first; + + StateParameters state_parameters_ = { .number_of_grid_cells_ = this->number_of_grid_cells_, + .number_of_species_ = number_of_species, + .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_ = diagonal_elements }; + + return Solver, State>( + BackwardEuler( + parameters, linear_solver, process_set, jacobian_diagonal_elements, this->reactions_), + state_parameters_, + this->number_of_grid_cells_, + number_of_species, + this->reactions_.size()); + } + +} // namespace micm \ No newline at end of file diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index fde718da7..4ff66d41d 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -35,25 +35,25 @@ namespace micm std::set> nonzero_jacobian_elements_{}; }; - template class MatrixPolicy = Matrix, template class SparseMatrixPolicy = StandardSparseMatrix> + template struct State { /// @brief The concentration of chemicals, varies through time - MatrixPolicy variables_; + DenseMatrixPolicy variables_; /// @brief Rate paramters particular to user-defined rate constants, may vary in time - MatrixPolicy custom_rate_parameters_; + DenseMatrixPolicy custom_rate_parameters_; /// @brief The reaction rates, may vary in time - MatrixPolicy rate_constants_; + DenseMatrixPolicy rate_constants_; /// @brief Atmospheric conditions, varies in time std::vector conditions_; /// @brief The jacobian structure, varies for each solve - SparseMatrixPolicy jacobian_; + SparseMatrixPolicy jacobian_; /// @brief Immutable data required for the state std::map variable_map_; std::map custom_rate_parameter_map_; std::vector variable_names_{}; - SparseMatrixPolicy lower_matrix_; - SparseMatrixPolicy upper_matrix_; + SparseMatrixPolicy lower_matrix_; + SparseMatrixPolicy upper_matrix_; std::size_t state_size_; std::size_t number_of_grid_cells_; diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 164e48c3f..186b7d0f8 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -58,8 +58,8 @@ inline std::error_code make_error_code(MicmStateErrc e) namespace micm { - template class MatrixPolicy, template class SparseMatrixPolicy> - inline State::State() + template + inline State::State() : conditions_(), variables_(), custom_rate_parameters_(), @@ -68,8 +68,8 @@ namespace micm { } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline State::State(const StateParameters& parameters) + template + inline State::State(const StateParameters& parameters) : conditions_(parameters.number_of_grid_cells_), variables_(parameters.number_of_grid_cells_, parameters.variable_names_.size(), 0.0), custom_rate_parameters_(parameters.number_of_grid_cells_, parameters.custom_rate_parameter_labels_.size(), 0.0), @@ -93,15 +93,15 @@ namespace micm jacobian_ = BuildJacobian( parameters.nonzero_jacobian_elements_, parameters.number_of_grid_cells_, state_size_); - auto lu = LuDecomposition::GetLUMatrices(jacobian_, 1.0e-30); + auto lu = LuDecomposition::GetLUMatrices(jacobian_, 1.0e-30); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); lower_matrix_ = lower_matrix; upper_matrix_ = upper_matrix; } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::SetConcentrations( + template + inline void State::SetConcentrations( const std::unordered_map>& species_to_concentration) { const std::size_t num_grid_cells = conditions_.size(); @@ -109,8 +109,8 @@ namespace micm SetConcentration({ pair.first }, pair.second); } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::SetConcentration(const Species& species, double concentration) + template + inline void State::SetConcentration(const Species& species, double concentration) { auto var = variable_map_.find(species.name_); if (var == variable_map_.end()) @@ -120,8 +120,8 @@ namespace micm variables_[0][variable_map_[species.name_]] = concentration; } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::SetConcentration( + template + inline void State::SetConcentration( const Species& species, const std::vector& concentration) { @@ -135,8 +135,8 @@ namespace micm variables_[i][i_species] = concentration[i]; } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::UnsafelySetCustomRateParameters( + template + inline void State::UnsafelySetCustomRateParameters( const std::vector>& parameters) { if (parameters.size() != variables_.NumRows()) @@ -152,16 +152,16 @@ namespace micm } } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::SetCustomRateParameters( + template + inline void State::SetCustomRateParameters( const std::unordered_map>& parameters) { for (auto& pair : parameters) SetCustomRateParameter(pair.first, pair.second); } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::SetCustomRateParameter(const std::string& label, double value) + template + inline void State::SetCustomRateParameter(const std::string& label, double value) { auto param = custom_rate_parameter_map_.find(label); if (param == custom_rate_parameter_map_.end()) @@ -172,8 +172,8 @@ namespace micm custom_rate_parameters_[0][param->second] = value; } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::SetCustomRateParameter( + template + inline void State::SetCustomRateParameter( const std::string& label, const std::vector& values) { @@ -187,8 +187,8 @@ namespace micm custom_rate_parameters_[i][param->second] = values[i]; } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::PrintHeader() + template + inline void State::PrintHeader() { auto largest_str_iter = std::max_element( variable_names_.begin(), variable_names_.end(), [](const auto& a, const auto& b) { return a.size() < b.size(); }); @@ -208,8 +208,8 @@ namespace micm std::cout << std::endl; } - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void State::PrintState(double time) + template + inline void State::PrintState(double time) { std::ios oldState(nullptr); oldState.copyfmt(std::cout); diff --git a/include/micm/util/cuda_dense_matrix.hpp b/include/micm/util/cuda_dense_matrix.hpp index 5b9a00ad5..41a5c517a 100644 --- a/include/micm/util/cuda_dense_matrix.hpp +++ b/include/micm/util/cuda_dense_matrix.hpp @@ -62,6 +62,11 @@ namespace micm template class CudaDenseMatrix : public VectorMatrix { + public: + // Diagonal markowitz reordering requires an int argument, make sure one is always accessible + using IntMatrix = CudaDenseMatrix; + using value_type = T; + private: /// @brief The device pointer (handle) to the allocated memory on the target device. CudaMatrixParam param_; diff --git a/include/micm/util/cuda_sparse_matrix.hpp b/include/micm/util/cuda_sparse_matrix.hpp index 5de9a534e..aa0a6bfd4 100644 --- a/include/micm/util/cuda_sparse_matrix.hpp +++ b/include/micm/util/cuda_sparse_matrix.hpp @@ -20,6 +20,11 @@ namespace micm template class CudaSparseMatrix : public SparseMatrix { + public: + // Diagonal markowitz reordering requires an int argument, make sure one is always accessible + using IntMatrix = CudaSparseMatrix; + using value_type = T; + private: CudaMatrixParam param_; diff --git a/include/micm/util/error.hpp b/include/micm/util/error.hpp index cdc0f4f73..145bc42e2 100644 --- a/include/micm/util/error.hpp +++ b/include/micm/util/error.hpp @@ -25,6 +25,7 @@ #define MICM_ERROR_CATEGORY_JIT "MICM JIT" #define MICM_JIT_ERROR_CODE_INVALID_MATRIX 1 #define MICM_JIT_ERROR_CODE_MISSING_JIT_FUNCTION 2 +#define MICM_JIT_ERROR_CODE_FAILED_TO_BUILD 3 #define MICM_ERROR_CATEGORY_PROCESS "MICM Process" #define MICM_PROCESS_ERROR_CODE_TOO_MANY_REACTANTS_FOR_SURFACE_REACTION 1 diff --git a/include/micm/util/jacobian.hpp b/include/micm/util/jacobian.hpp index 444953d4b..2cad50cb6 100644 --- a/include/micm/util/jacobian.hpp +++ b/include/micm/util/jacobian.hpp @@ -13,15 +13,15 @@ namespace micm { // annonymous namespace to hide jacobian builder - template class SparseMatrixPolicy> - SparseMatrixPolicy BuildJacobian( + template + SparseMatrixPolicy BuildJacobian( const std::set>& nonzero_jacobian_elements, std::size_t number_of_grid_cells, std::size_t state_size) { MICM_PROFILE_FUNCTION(); - auto builder = SparseMatrixPolicy::Create(state_size).SetNumberOfBlocks(number_of_grid_cells); + auto builder = SparseMatrixPolicy::Create(state_size).SetNumberOfBlocks(number_of_grid_cells); for (auto& elem : nonzero_jacobian_elements) builder = builder.WithElement(elem.first, elem.second); // Always include diagonal elements @@ -29,6 +29,6 @@ namespace micm { builder = builder.WithElement(i, i); } - return SparseMatrixPolicy(builder); + return SparseMatrixPolicy(builder); } } // namespace micm diff --git a/include/micm/util/matrix.hpp b/include/micm/util/matrix.hpp index 3d044e7fb..0b2b4c2ee 100644 --- a/include/micm/util/matrix.hpp +++ b/include/micm/util/matrix.hpp @@ -25,9 +25,15 @@ namespace micm }; /// @brief A 2D array class with contiguous memory - template + template class Matrix { + public: + // Diagonal markowitz reordering requires an int argument, make sure one is always accessible + using IntMatrix = Matrix; + using value_type = T; + + private: std::vector data_; std::size_t x_dim_; std::size_t y_dim_; @@ -274,4 +280,6 @@ namespace micm } }; + using StandardDenseMatrix = Matrix; + } // namespace micm diff --git a/include/micm/util/sparse_matrix.hpp b/include/micm/util/sparse_matrix.hpp index 5ed2a6285..92470fb9e 100644 --- a/include/micm/util/sparse_matrix.hpp +++ b/include/micm/util/sparse_matrix.hpp @@ -32,8 +32,7 @@ namespace micm template class SparseMatrix; - template - using StandardSparseMatrix = SparseMatrix; + using StandardSparseMatrix = SparseMatrix; /// @brief A sparse block-diagonal 2D matrix class with contiguous memory /// @@ -43,9 +42,14 @@ namespace micm /// /// The template parameters are the type of the matrix elements and a class that /// defines the sizing and ordering of the data elements - template + template class SparseMatrix : public OrderingPolicy { + public: + // Diagonal markowitz reordering requires an int argument, make sure one is always accessible + using IntMatrix = SparseMatrix; + using value_type = T; + protected: std::size_t number_of_blocks_; // Number of block sub-matrices in the overall matrix std::vector row_ids_; // Row indices of each non-zero element in a block diff --git a/include/micm/util/sparse_matrix_vector_ordering.hpp b/include/micm/util/sparse_matrix_vector_ordering.hpp index 50c09b561..58b59a0e9 100644 --- a/include/micm/util/sparse_matrix_vector_ordering.hpp +++ b/include/micm/util/sparse_matrix_vector_ordering.hpp @@ -69,7 +69,6 @@ namespace micm }; // Default vectorized SparseMatrix - template - using VectorSparseMatrix = SparseMatrix>; + using DefaultVectorSparseMatrix = SparseMatrix>; } // namespace micm diff --git a/include/micm/util/vector_matrix.hpp b/include/micm/util/vector_matrix.hpp index 989868edb..12a368024 100644 --- a/include/micm/util/vector_matrix.hpp +++ b/include/micm/util/vector_matrix.hpp @@ -31,6 +31,12 @@ namespace micm template class VectorMatrix { + public: + // Diagonal markowitz reordering requires an int argument, make sure one is always accessible + using IntMatrix = VectorMatrix; + using value_type = T; + + private: protected: std::vector data_; std::size_t x_dim_; // number of rows diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3d2c33410..a6eb78abf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -99,6 +99,7 @@ if(MICM_ENABLE_CUDA) PUBLIC cudart cublas culibos ) + target_include_directories(micm_cuda PUBLIC /usr/local/cuda/include) endif() if (MICM_ENABLE_PROFILE) diff --git a/test/integration/analytical_jit_rosenbrock.cpp b/test/integration/analytical_jit_rosenbrock.cpp index 9d40b986a..65ae4452e 100644 --- a/test/integration/analytical_jit_rosenbrock.cpp +++ b/test/integration/analytical_jit_rosenbrock.cpp @@ -8,8 +8,7 @@ template using DefaultVectorMatrix = micm::VectorMatrix; -template -using DefaultSparseVectorMatrix = micm::SparseMatrix>; +using DefaultSparseVectorMatrix = micm::SparseMatrix>; using DefaultJitRosenbrockSolver = micm::JitRosenbrockSolver< DefaultVectorMatrix, @@ -19,238 +18,70 @@ using DefaultJitRosenbrockSolver = micm::JitRosenbrockSolver< TEST(AnalyticalExamplesJitRosenbrock, Troe) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_troe( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_troe(); } TEST(AnalyticalExamplesJitRosenbrock, TroeSuperStiffButAnalytical) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_stiff_troe( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_troe(); } TEST(AnalyticalExamplesJitRosenbrock, Photolysis) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_photolysis( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_photolysis(); } TEST(AnalyticalExamplesJitRosenbrock, PhotolysisSuperStiffButAnalytical) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_stiff_photolysis( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_photolysis(); } TEST(AnalyticalExamplesJitRosenbrock, TernaryChemicalActivation) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_ternary_chemical_activation( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_ternary_chemical_activation(); } TEST(AnalyticalExamplesJitRosenbrock, TernaryChemicalActivationSuperStiffButAnalytical) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_stiff_ternary_chemical_activation( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_ternary_chemical_activation(); } TEST(AnalyticalExamplesJitRosenbrock, Tunneling) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_tunneling( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_tunneling(); } TEST(AnalyticalExamplesJitRosenbrock, TunnelingSuperStiffButAnalytical) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_stiff_tunneling( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_tunneling(); } TEST(AnalyticalExamplesJitRosenbrock, Arrhenius) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_arrhenius( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_arrhenius(); } TEST(AnalyticalExamplesJitRosenbrock, ArrheniusSuperStiffButAnalytical) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_stiff_arrhenius( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_arrhenius(); } TEST(AnalyticalExamplesJitRosenbrock, Branched) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_branched( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_branched(); } TEST(AnalyticalExamplesJitRosenbrock, BranchedSuperStiffButAnalytical) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_stiff_branched( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_branched(); } TEST(AnalyticalExamplesJitRosenbrock, Robertson) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_robertson( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_robertson(); } TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - test_analytical_surface_rxn( - [&](const micm::System& s, const std::vector& p) -> DefaultJitRosenbrockSolver - { - return DefaultJitRosenbrockSolver{ - jit.get(), s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_surface_rxn(); } diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp index f0de34075..a44a8202b 100644 --- a/test/integration/analytical_policy.hpp +++ b/test/integration/analytical_policy.hpp @@ -71,12 +71,11 @@ void writeCSV( using yields = std::pair; -template -using SparseMatrixTest = micm::SparseMatrix; +using SparseMatrixTest = micm::SparseMatrix; template void test_analytical_troe( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A -> B, k1 @@ -111,9 +110,8 @@ void test_analytical_troe( .SetPhase(gas_phase); auto processes = std::vector{ r1, r2 }; - OdeSolverPolicy solver = create_solver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), processes); + OdeSolverPolicy solver = OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), processes, parameters); - micm::BackwardEuler be; auto be_state = solver.GetState(); double temperature = 272.5; @@ -162,18 +160,10 @@ void test_analytical_troe( times.push_back(time_step); // Model results auto result = solver.Solve(time_step, state); - if constexpr (std::is_same_v>) - { - auto linear_solver = solver.linear_solver_; - auto process_set = solver.process_set_; - be.Solve( - time_step, be_state, linear_solver, process_set, processes, solver.state_parameters_.jacobian_diagonal_elements_); - } EXPECT_EQ(result.state_, (micm::SolverState::Converged)); EXPECT_NEAR(k1, state.rate_constants_.AsVector()[0], 1e-8); EXPECT_NEAR(k2, state.rate_constants_.AsVector()[1], 1e-8); model_concentrations[i_time] = result.result_.AsVector(); - be_model_concentrations[i_time] = be_state.variables_[0]; state.variables_[0] = result.result_.AsVector(); // Analytical results @@ -205,22 +195,12 @@ void test_analytical_troe( << "Arrays differ at index (" << i << ", " << 1 << ")"; EXPECT_NEAR(model_concentrations[i][_c], analytical_concentrations[i][2], 1e-8) << "Arrays differ at index (" << i << ", " << 2 << ")"; - - if constexpr (std::is_same_v>) - { - EXPECT_NEAR(be_model_concentrations[i][_a], analytical_concentrations[i][0], 1e-4) - << "Arrays differ at index (" << i << ", " << 0 << ")"; - EXPECT_NEAR(be_model_concentrations[i][_b], analytical_concentrations[i][1], 1e-4) - << "Arrays differ at index (" << i << ", " << 1 << ")"; - EXPECT_NEAR(be_model_concentrations[i][_c], analytical_concentrations[i][2], 1e-4) - << "Arrays differ at index (" << i << ", " << 2 << ")"; - } } } template void test_analytical_stiff_troe( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A1 -> B, k1 @@ -275,8 +255,8 @@ void test_analytical_stiff_troe( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 0.9 * 4.0e10 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -358,7 +338,7 @@ void test_analytical_stiff_troe( template void test_analytical_photolysis( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A -> B, k1 @@ -386,7 +366,7 @@ void test_analytical_photolysis( .SetPhase(gas_phase); OdeSolverPolicy solver = - create_solver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }); + OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -464,7 +444,7 @@ void test_analytical_photolysis( template void test_analytical_stiff_photolysis( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A1 -> B, k1 @@ -512,8 +492,8 @@ void test_analytical_stiff_photolysis( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 0.9 * 4.0e10 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -594,7 +574,7 @@ void test_analytical_stiff_photolysis( template void test_analytical_ternary_chemical_activation( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A -> B, k1 @@ -630,7 +610,7 @@ void test_analytical_ternary_chemical_activation( .SetPhase(gas_phase); OdeSolverPolicy solver = - create_solver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }); + OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -711,7 +691,7 @@ void test_analytical_ternary_chemical_activation( template void test_analytical_stiff_ternary_chemical_activation( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A1 -> B, k1 @@ -766,8 +746,8 @@ void test_analytical_stiff_ternary_chemical_activation( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 0.9 * 4.0e10 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -849,7 +829,7 @@ void test_analytical_stiff_ternary_chemical_activation( template void test_analytical_tunneling( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A -> B, k1 @@ -878,7 +858,7 @@ void test_analytical_tunneling( .SetPhase(gas_phase); OdeSolverPolicy solver = - create_solver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }); + OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -953,7 +933,7 @@ void test_analytical_tunneling( template void test_analytical_stiff_tunneling( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A1 -> B, k1 @@ -1001,8 +981,8 @@ void test_analytical_stiff_tunneling( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 0.9 * 4.0e10 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1081,7 +1061,7 @@ void test_analytical_stiff_tunneling( template void test_analytical_arrhenius( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A -> B, k1 @@ -1109,7 +1089,7 @@ void test_analytical_arrhenius( .SetPhase(gas_phase); OdeSolverPolicy solver = - create_solver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }); + OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1184,7 +1164,7 @@ void test_analytical_arrhenius( template void test_analytical_stiff_arrhenius( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A1 -> B, k1 @@ -1233,8 +1213,8 @@ void test_analytical_stiff_arrhenius( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 0.9 * 4.0e10 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1313,7 +1293,7 @@ void test_analytical_stiff_arrhenius( template void test_analytical_branched( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A -> B, k1 @@ -1351,7 +1331,7 @@ void test_analytical_branched( .SetPhase(gas_phase); OdeSolverPolicy solver = - create_solver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }); + OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1441,7 +1421,7 @@ void test_analytical_branched( template void test_analytical_stiff_branched( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A1 -> B, k1 @@ -1504,8 +1484,8 @@ void test_analytical_stiff_branched( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 0.9 * 4.0e10 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1599,7 +1579,7 @@ void test_analytical_stiff_branched( template void test_analytical_robertson( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { /* * A -> B, k1 = 0.04 @@ -1638,8 +1618,8 @@ void test_analytical_robertson( .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r3" })) .SetPhase(gas_phase); - OdeSolverPolicy solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3 }); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3 }, parameters); double temperature = 272.5; double pressure = 101253.3; diff --git a/test/integration/analytical_rosenbrock.cpp b/test/integration/analytical_rosenbrock.cpp index f3b2843d3..5383b649a 100644 --- a/test/integration/analytical_rosenbrock.cpp +++ b/test/integration/analytical_rosenbrock.cpp @@ -9,175 +9,76 @@ #include -template -using SparseMatrixTest = micm::SparseMatrix; +using SparseMatrixTest = micm::SparseMatrix; TEST(AnalyticalExamples, Troe) { - test_analytical_troe>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_troe>(); } TEST(AnalyticalExamples, TroeSuperStiffButAnalytical) { - test_analytical_stiff_troe>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_troe>(); } TEST(AnalyticalExamples, Photolysis) { - test_analytical_photolysis>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_photolysis>(); } TEST(AnalyticalExamples, PhotolysisSuperStiffButAnalytical) { - test_analytical_stiff_photolysis>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_photolysis>(); } TEST(AnalyticalExamples, TernaryChemicalActivation) { - test_analytical_ternary_chemical_activation>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_ternary_chemical_activation>(); } TEST(AnalyticalExamples, TernaryChemicalActivationSuperStiffButAnalytical) { - test_analytical_stiff_ternary_chemical_activation>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_ternary_chemical_activation>(); } TEST(AnalyticalExamples, Tunneling) { - test_analytical_tunneling>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_tunneling>(); } TEST(AnalyticalExamples, TunnelingSuperStiffButAnalytical) { - test_analytical_stiff_tunneling>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_tunneling>(); } TEST(AnalyticalExamples, Arrhenius) { - test_analytical_arrhenius>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_arrhenius>(); } TEST(AnalyticalExamples, ArrheniusSuperStiffButAnalytical) { - test_analytical_stiff_arrhenius>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_arrhenius>(); } TEST(AnalyticalExamples, Branched) { - test_analytical_branched>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_branched>(); } TEST(AnalyticalExamples, BranchedSuperStiffButAnalytical) { - test_analytical_stiff_branched>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_stiff_branched>(); } TEST(AnalyticalExamples, Robertson) { - test_analytical_robertson>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_robertson>(); } TEST(AnalyticalExamples, SurfaceRxn) { - test_analytical_surface_rxn>( - [](const micm::System& s, - const std::vector& p) -> micm::RosenbrockSolver - { - return micm::RosenbrockSolver{ - s, p, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters() - }; - }); + test_analytical_surface_rxn>(); } TEST(AnalyticalExamples, Oregonator) diff --git a/test/integration/analytical_surface_rxn_policy.hpp b/test/integration/analytical_surface_rxn_policy.hpp index 452cefb4e..6547496a0 100644 --- a/test/integration/analytical_surface_rxn_policy.hpp +++ b/test/integration/analytical_surface_rxn_policy.hpp @@ -9,7 +9,7 @@ template void test_analytical_surface_rxn( - const std::function&)> create_solver) + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { // parameters, from CAMP/test/unit_rxn_data/test_rxn_surface.F90 const double mode_GMD = 1.0e-6; // mode geometric mean diameter [m] @@ -69,7 +69,7 @@ void test_analytical_surface_rxn( // }; // Solver - OdeSolverPolicy solver = create_solver(chemical_system, reactions); + OdeSolverPolicy solver = OdeSolverPolicy(chemical_system, reactions, parameters); // State micm::State state = solver.GetState(); diff --git a/test/integration/chapman.cpp b/test/integration/chapman.cpp index 6b1df9329..6b6c24df9 100644 --- a/test/integration/chapman.cpp +++ b/test/integration/chapman.cpp @@ -14,8 +14,7 @@ using yields = std::pair; -template -using SparseMatrixTest = micm::SparseMatrix; +using SparseMatrixTest = micm::SparseMatrix; #ifdef USE_JSON # include diff --git a/test/integration/cmake/test_micm.cpp b/test/integration/cmake/test_micm.cpp index ccd02e835..0310bbc00 100644 --- a/test/integration/cmake/test_micm.cpp +++ b/test/integration/cmake/test_micm.cpp @@ -11,8 +11,8 @@ void print_header() << "," << std::setw(10) << "C" << std::endl; } -template class T> -void print_state(double time, State& state) +template +void print_state(double time, StatePolicy& state) { std::ios oldState(nullptr); oldState.copyfmt(std::cout); @@ -42,7 +42,7 @@ int main() auto reactions = solver_params.processes_; RosenbrockSolver<> solver{ chemical_system, reactions, params }; - State state = solver.GetState(); + auto state = solver.GetState(); // mol m-3 state.variables_[0] = { 1, 0, 0 }; diff --git a/test/integration/e5.hpp b/test/integration/e5.hpp index 54caa9a7d..db0a771d6 100644 --- a/test/integration/e5.hpp +++ b/test/integration/e5.hpp @@ -1,11 +1,12 @@ #pragma once #include +#include template< template class MatrixPolicy = micm::Matrix, - template class SparseMatrixPolicy = micm::SparseMatrix, - class LinearSolverPolicy = micm::LinearSolver> + class SparseMatrixPolicy = micm::StandardSparseMatrix, + class LinearSolverPolicy = micm::LinearSolver> class E5 : public micm::RosenbrockSolver { std::set> nonzero_jacobian_elements_; @@ -22,7 +23,7 @@ class E5 : public micm::RosenbrockSolverprocesses_ = processes; this->parameters_ = parameters; - auto builder = SparseMatrixPolicy::Create(4).SetNumberOfBlocks(1).InitialValue(0.0); + auto builder = SparseMatrixPolicy::Create(4).SetNumberOfBlocks(1).InitialValue(0.0); for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) @@ -31,7 +32,7 @@ class E5 : public micm::RosenbrockSolver jacobian = SparseMatrixPolicy(builder); + SparseMatrixPolicy jacobian = SparseMatrixPolicy(builder); std::vector jacobian_diagonal_elements; for (std::size_t i = 0; i < jacobian.NumRows(); ++i) @@ -60,9 +61,9 @@ class E5 : public micm::RosenbrockSolver GetState() const override + micm::State, SparseMatrixPolicy> GetState() const override { - auto state = micm::State{ this->state_parameters_ }; + auto state = micm::State, SparseMatrixPolicy>{ this->state_parameters_ }; state.jacobian_ = micm::BuildJacobian( nonzero_jacobian_elements_, @@ -108,7 +109,7 @@ class E5 : public micm::RosenbrockSolver& rate_constants, const MatrixPolicy& number_densities, - SparseMatrixPolicy& jacobian) override + SparseMatrixPolicy& jacobian) override { std::fill(jacobian.AsVector().begin(), jacobian.AsVector().end(), 0.0); auto data = number_densities.AsVector(); diff --git a/test/integration/hires.hpp b/test/integration/hires.hpp index 73f376143..ac1ccd83f 100644 --- a/test/integration/hires.hpp +++ b/test/integration/hires.hpp @@ -4,8 +4,8 @@ template< template class MatrixPolicy = micm::Matrix, - template class SparseMatrixPolicy = micm::SparseMatrix, - class LinearSolverPolicy = micm::LinearSolver> + class SparseMatrixPolicy = micm::StandardSparseMatrix, + class LinearSolverPolicy = micm::LinearSolver> class HIRES : public micm::RosenbrockSolver { std::set> nonzero_jacobian_elements_; @@ -23,7 +23,7 @@ class HIRES : public micm::RosenbrockSolverprocesses_ = processes; this->parameters_ = parameters; - auto builder = SparseMatrixPolicy::Create(8).SetNumberOfBlocks(1).InitialValue(0.0); + auto builder = SparseMatrixPolicy::Create(8).SetNumberOfBlocks(1).InitialValue(0.0); for (int i = 0; i < 8; ++i) { for (int j = 0; j < 8; ++j) @@ -32,7 +32,7 @@ class HIRES : public micm::RosenbrockSolver jacobian = SparseMatrixPolicy(builder); + SparseMatrixPolicy jacobian = SparseMatrixPolicy(builder); std::vector jacobian_diagonal_elements; for (std::size_t i = 0; i < jacobian.NumRows(); ++i) @@ -61,9 +61,9 @@ class HIRES : public micm::RosenbrockSolver GetState() const override + micm::State, SparseMatrixPolicy> GetState() const override { - auto state = micm::State{ this->state_parameters_ }; + auto state = micm::State, SparseMatrixPolicy>{ this->state_parameters_ }; state.jacobian_ = micm::BuildJacobian( nonzero_jacobian_elements_, @@ -109,7 +109,7 @@ class HIRES : public micm::RosenbrockSolver& rate_constants, const MatrixPolicy& number_densities, - SparseMatrixPolicy& jacobian) override + SparseMatrixPolicy& jacobian) override { std::fill(jacobian.AsVector().begin(), jacobian.AsVector().end(), 0.0); auto data = number_densities.AsVector(); diff --git a/test/integration/integrated_reaction_rates.cpp b/test/integration/integrated_reaction_rates.cpp index 698d3ba6f..06c541012 100644 --- a/test/integration/integrated_reaction_rates.cpp +++ b/test/integration/integrated_reaction_rates.cpp @@ -14,8 +14,7 @@ using yields = std::pair; -template -using SparseMatrixTest = micm::SparseMatrix; +using SparseMatrixTest = micm::SparseMatrix; TEST(ChapmanIntegration, CanBuildChapmanSystem) { diff --git a/test/integration/jit_terminator.cpp b/test/integration/jit_terminator.cpp index c8c42ec5e..254e6c5fa 100644 --- a/test/integration/jit_terminator.cpp +++ b/test/integration/jit_terminator.cpp @@ -9,39 +9,20 @@ #include -template class MatrixPolicy, template class SparseMatrixPolicy> +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< - MatrixPolicy, micm::JitRosenbrockSolver< MatrixPolicy, SparseMatrixPolicy, micm::JitLinearSolver, micm::JitProcessSet>>( - [&](const micm::System& s, const std::vector& p) - -> micm::JitRosenbrockSolver< - MatrixPolicy, - SparseMatrixPolicy, - micm::JitLinearSolver, - micm::JitProcessSet> - { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - 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; - return micm::JitRosenbrockSolver< - MatrixPolicy, - SparseMatrixPolicy, - micm::JitLinearSolver, - micm::JitProcessSet>{ jit.get(), s, p, solver_params }; - }, - number_of_grid_cells); + number_of_grid_cells, solver_params); } template @@ -53,14 +34,10 @@ using Group3VectorMatrix = micm::VectorMatrix; template using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; TEST(JitRosenbrockSolver, Terminator) { diff --git a/test/integration/oregonator.hpp b/test/integration/oregonator.hpp index d8117973f..7ccc857c6 100644 --- a/test/integration/oregonator.hpp +++ b/test/integration/oregonator.hpp @@ -4,8 +4,8 @@ template< template class MatrixPolicy = micm::Matrix, - template class SparseMatrixPolicy = micm::SparseMatrix, - class LinearSolverPolicy = micm::LinearSolver> + class SparseMatrixPolicy = micm::StandardSparseMatrix, + class LinearSolverPolicy = micm::LinearSolver> class Oregonator : public micm::RosenbrockSolver { std::set> nonzero_jacobian_elements_; @@ -24,7 +24,7 @@ class Oregonator : public micm::RosenbrockSolverparameters_ = parameters; this->parameters_.absolute_tolerance_ = std::vector(3, 1.0e-3); - auto builder = SparseMatrixPolicy::Create(3).SetNumberOfBlocks(1).InitialValue(0.0); + auto builder = SparseMatrixPolicy::Create(3).SetNumberOfBlocks(1).InitialValue(0.0); for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) @@ -33,7 +33,7 @@ class Oregonator : public micm::RosenbrockSolver jacobian = SparseMatrixPolicy(builder); + SparseMatrixPolicy jacobian = SparseMatrixPolicy(builder); std::vector jacobian_diagonal_elements; for (std::size_t i = 0; i < jacobian.NumRows(); ++i) @@ -62,9 +62,9 @@ class Oregonator : public micm::RosenbrockSolver GetState() const override + micm::State, SparseMatrixPolicy> GetState() const override { - auto state = micm::State{ this->state_parameters_ }; + auto state = micm::State, SparseMatrixPolicy>{ this->state_parameters_ }; state.jacobian_ = micm::BuildJacobian( nonzero_jacobian_elements_, @@ -105,7 +105,7 @@ class Oregonator : public micm::RosenbrockSolver& rate_constants, const MatrixPolicy& number_densities, - SparseMatrixPolicy& jacobian) override + SparseMatrixPolicy& jacobian) override { auto data = number_densities.AsVector(); diff --git a/test/integration/terminator.cpp b/test/integration/terminator.cpp index fc3256aa3..0027e3463 100644 --- a/test/integration/terminator.cpp +++ b/test/integration/terminator.cpp @@ -10,41 +10,27 @@ #include -template -using SparseMatrixTest = micm::SparseMatrix; +using SparseMatrixTest = micm::SparseMatrix; -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> void RunTerminatorTest(std::size_t number_of_grid_cells) { - TestTerminator>( - [&](const micm::System& s, const std::vector& p) - -> micm::RosenbrockSolver - { - 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; - return micm::RosenbrockSolver{ s, p, solver_params }; - }, - number_of_grid_cells); - TestTerminator>( - [&](const micm::System& s, const std::vector& p) - -> micm::RosenbrockSolver - { - 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; - solver_params.check_singularity_ = true; - return micm::RosenbrockSolver{ s, p, solver_params }; - }, - number_of_grid_cells); + 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>(number_of_grid_cells, solver_params); + + solver_params.check_singularity_ = true; + TestTerminator< micm::RosenbrockSolver>(number_of_grid_cells, solver_params); } TEST(RosenbrockSolver, Terminator) { - RunTerminatorTest>(2); - RunTerminatorTest>(2); - RunTerminatorTest>(3); - RunTerminatorTest>(4); + RunTerminatorTest>(2); + RunTerminatorTest>(2); + RunTerminatorTest>(3); + RunTerminatorTest>(4); } template @@ -56,19 +42,15 @@ using Group3VectorMatrix = micm::VectorMatrix; template using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; TEST(RosenbrockSolver, VectorTerminator) { - RunTerminatorTest>(1); - RunTerminatorTest>(4); - RunTerminatorTest>(3); - RunTerminatorTest>(2); + RunTerminatorTest>(1); + RunTerminatorTest>(4); + RunTerminatorTest>(3); + RunTerminatorTest>(2); } \ No newline at end of file diff --git a/test/integration/terminator.hpp b/test/integration/terminator.hpp index 0dc58410e..7a14babfb 100644 --- a/test/integration/terminator.hpp +++ b/test/integration/terminator.hpp @@ -20,10 +20,10 @@ /// /// More details including analytical solution can be found here: /// https://github.com/ESCOMP/CAM/blob/8cd44c50fe107c0b93ccd48b61eaa3d10a5b4e2f/src/chemistry/pp_terminator/chemistry.F90#L1-L434 -template class MatrixPolicy, class OdeSolverPolicy> +template void TestTerminator( - const std::function&)> create_solver, - std::size_t number_of_grid_cells) + std::size_t number_of_grid_cells, + const micm::RosenbrockSolverParameters parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) { auto cl2 = micm::Species("Cl2"); auto cl = micm::Species("Cl"); @@ -45,8 +45,8 @@ void TestTerminator( .SetPhase(gas_phase) .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = k2 })); - auto solver = create_solver( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ toy_r1, toy_r2 }); + auto solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ toy_r1, toy_r2 }, parameters); auto state = solver.GetState(); auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); diff --git a/test/regression/RosenbrockChapman/jit_util.hpp b/test/regression/RosenbrockChapman/jit_util.hpp index 0e0fd3369..fe9e1caba 100644 --- a/test/regression/RosenbrockChapman/jit_util.hpp +++ b/test/regression/RosenbrockChapman/jit_util.hpp @@ -5,7 +5,6 @@ template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -18,20 +17,12 @@ getTwoStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - return micm::JitRosenbrockSolver( - jit.get(), micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -44,20 +35,12 @@ getThreeStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - return micm::JitRosenbrockSolver( - jit.get(), micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -70,20 +53,12 @@ getFourStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::FourStageRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - return micm::JitRosenbrockSolver( - jit.get(), micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -96,20 +71,12 @@ getFourStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - return micm::JitRosenbrockSolver( - jit.get(), micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } template< template class MatrixPolicy, - template class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> @@ -122,12 +89,5 @@ getSixStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - return micm::JitRosenbrockSolver( - jit.get(), micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver( micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_dforce_dy.cpp b/test/regression/RosenbrockChapman/regression_test_dforce_dy.cpp index 06cec436c..aeebd280f 100644 --- a/test/regression/RosenbrockChapman/regression_test_dforce_dy.cpp +++ b/test/regression/RosenbrockChapman/regression_test_dforce_dy.cpp @@ -6,7 +6,7 @@ TEST(RegressionRosenbrock, Jacobian) { - auto solver = getThreeStageMultiCellChapmanSolver>(3); + auto solver = getThreeStageMultiCellChapmanSolver>(3); testJacobian<>(solver); } @@ -16,28 +16,28 @@ TEST(RegressionRosenbrock, VectorJacobian) auto solver = getThreeStageMultiCellChapmanSolver< Group1VectorMatrix, Group1SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testJacobian<>(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group2VectorMatrix, Group2SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testJacobian<>(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group3VectorMatrix, Group3SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testJacobian<>(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group4VectorMatrix, Group4SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testJacobian<>(solver); } } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_dforce_dy_policy.hpp b/test/regression/RosenbrockChapman/regression_test_dforce_dy_policy.hpp index 44c7574f4..b9b5571ed 100644 --- a/test/regression/RosenbrockChapman/regression_test_dforce_dy_policy.hpp +++ b/test/regression/RosenbrockChapman/regression_test_dforce_dy_policy.hpp @@ -51,8 +51,7 @@ void testJacobian(OdeSolverPolicy& solver) template using DenseMatrix = micm::Matrix; -template -using SparseMatrix = micm::SparseMatrix; +using SparseMatrix = micm::SparseMatrix; template using Group1VectorMatrix = micm::VectorMatrix; @@ -63,11 +62,7 @@ using Group3VectorMatrix = micm::VectorMatrix; template using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; \ No newline at end of file +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_p_force.cpp b/test/regression/RosenbrockChapman/regression_test_p_force.cpp index 4be94cbb4..280a0abf5 100644 --- a/test/regression/RosenbrockChapman/regression_test_p_force.cpp +++ b/test/regression/RosenbrockChapman/regression_test_p_force.cpp @@ -6,7 +6,7 @@ TEST(RegressionRosenbrock, RateConstants) { - auto solver = getThreeStageMultiCellChapmanSolver>(3); + auto solver = getThreeStageMultiCellChapmanSolver>(3); testRateConstants<>(solver); } @@ -16,35 +16,35 @@ TEST(RegressionRosenbrock, VectorRateConstants) auto solver = getThreeStageMultiCellChapmanSolver< Group1VectorMatrix, Group1SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testRateConstants<>(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group2VectorMatrix, Group2SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testRateConstants<>(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group3VectorMatrix, Group3SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testRateConstants<>(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group4VectorMatrix, Group4SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testRateConstants<>(solver); } } TEST(RegressionRosenbrock, Forcing) { - auto solver = getThreeStageMultiCellChapmanSolver>(3); + auto solver = getThreeStageMultiCellChapmanSolver>(3); testForcing(solver); } @@ -54,28 +54,28 @@ TEST(RegressionRosenbrock, VectorForcing) auto solver = getThreeStageMultiCellChapmanSolver< Group1VectorMatrix, Group1SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testForcing(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group2VectorMatrix, Group2SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testForcing(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group3VectorMatrix, Group3SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testForcing(solver); } { auto solver = getThreeStageMultiCellChapmanSolver< Group4VectorMatrix, Group4SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testForcing(solver); } } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_p_force_policy.hpp b/test/regression/RosenbrockChapman/regression_test_p_force_policy.hpp index e0160ede0..ae498427e 100644 --- a/test/regression/RosenbrockChapman/regression_test_p_force_policy.hpp +++ b/test/regression/RosenbrockChapman/regression_test_p_force_policy.hpp @@ -85,8 +85,7 @@ void testForcing(OdeSolverPolicy& solver) template using DenseMatrix = micm::Matrix; -template -using SparseMatrix = micm::SparseMatrix; +using SparseMatrix = micm::SparseMatrix; template using Group1VectorMatrix = micm::VectorMatrix; @@ -97,11 +96,7 @@ using Group3VectorMatrix = micm::VectorMatrix; template using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; diff --git a/test/regression/RosenbrockChapman/regression_test_solve.cpp b/test/regression/RosenbrockChapman/regression_test_solve.cpp index 5c68f1885..7e1e4f0f1 100644 --- a/test/regression/RosenbrockChapman/regression_test_solve.cpp +++ b/test/regression/RosenbrockChapman/regression_test_solve.cpp @@ -6,31 +6,31 @@ TEST(RegressionRosenbrock, TwoStageSolve) { - auto solver = getTwoStageMultiCellChapmanSolver>(3); + auto solver = getTwoStageMultiCellChapmanSolver>(3); testSolve(solver, 1.0e-2); } TEST(RegressionRosenbrock, ThreeStageSolve) { - auto solver = getThreeStageMultiCellChapmanSolver>(3); + auto solver = getThreeStageMultiCellChapmanSolver>(3); testSolve(solver); } TEST(RegressionRosenbrock, FourStageSolve) { - auto solver = getFourStageMultiCellChapmanSolver>(3); + auto solver = getFourStageMultiCellChapmanSolver>(3); testSolve(solver, 1.0e-4); } TEST(RegressionRosenbrock, FourStageDASolve) { - auto solver = getFourStageDAMultiCellChapmanSolver>(3); + auto solver = getFourStageDAMultiCellChapmanSolver>(3); testSolve(solver, 1.0e-4); } TEST(RegressionRosenbrock, SixStageDASolve) { - auto solver = getSixStageDAMultiCellChapmanSolver>(3); + auto solver = getSixStageDAMultiCellChapmanSolver>(3); testSolve(solver, 1.0e-4); } @@ -39,24 +39,24 @@ TEST(RegressionRosenbrock, VectorSolve) auto solver1 = getThreeStageMultiCellChapmanSolver< Group1VectorMatrix, Group1SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testSolve(solver1); auto solver2 = getThreeStageMultiCellChapmanSolver< Group2VectorMatrix, Group2SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testSolve(solver2); auto solver3 = getThreeStageMultiCellChapmanSolver< Group3VectorMatrix, Group3SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testSolve(solver3); auto solver4 = getThreeStageMultiCellChapmanSolver< Group4VectorMatrix, Group4SparseVectorMatrix, - micm::LinearSolver>(3); + micm::LinearSolver>(3); testSolve(solver4); } \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/regression_test_solve_policy.hpp b/test/regression/RosenbrockChapman/regression_test_solve_policy.hpp index 1db2338b4..c06f5fe04 100644 --- a/test/regression/RosenbrockChapman/regression_test_solve_policy.hpp +++ b/test/regression/RosenbrockChapman/regression_test_solve_policy.hpp @@ -66,8 +66,8 @@ void testSolve(OdeSolverPolicy& solver, double relative_tolerance = 1.0e-8) template using DenseMatrix = micm::Matrix; -template -using SparseMatrix = micm::SparseMatrix; + +using SparseMatrix = micm::SparseMatrix; template using Group1VectorMatrix = micm::VectorMatrix; @@ -78,11 +78,7 @@ using Group3VectorMatrix = micm::VectorMatrix; template using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; \ No newline at end of file +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/util.hpp b/test/regression/RosenbrockChapman/util.hpp index 3d38456dd..72926e3fb 100644 --- a/test/regression/RosenbrockChapman/util.hpp +++ b/test/regression/RosenbrockChapman/util.hpp @@ -83,7 +83,7 @@ std::vector createProcesses(const micm::Phase& gas_phase) return { photo_1, photo_2, photo_3, r1, r2, r3, r4 }; } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> micm::RosenbrockSolver getTwoStageMultiCellChapmanSolver( const size_t number_of_grid_cells) { @@ -97,7 +97,7 @@ micm::RosenbrockSolver get micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> micm::RosenbrockSolver getThreeStageMultiCellChapmanSolver( const size_t number_of_grid_cells) { @@ -111,7 +111,7 @@ micm::RosenbrockSolver get micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> micm::RosenbrockSolver getFourStageMultiCellChapmanSolver( const size_t number_of_grid_cells) { @@ -125,7 +125,7 @@ micm::RosenbrockSolver get micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> micm::RosenbrockSolver getFourStageDAMultiCellChapmanSolver( const size_t number_of_grid_cells) { @@ -139,7 +139,7 @@ micm::RosenbrockSolver get micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> micm::RosenbrockSolver getSixStageDAMultiCellChapmanSolver( const size_t number_of_grid_cells) { diff --git a/test/tutorial/test_jit_tutorial.cpp b/test/tutorial/test_jit_tutorial.cpp index 4ea954804..fcee467ad 100644 --- a/test/tutorial/test_jit_tutorial.cpp +++ b/test/tutorial/test_jit_tutorial.cpp @@ -13,8 +13,7 @@ constexpr size_t n_grid_cells = 1; // partial template specializations template using GroupVectorMatrix = micm::VectorMatrix; -template -using GroupSparseVectorMatrix = micm::SparseMatrix>; +using GroupSparseVectorMatrix = micm::SparseMatrix>; auto run_solver(auto& solver) { @@ -94,15 +93,13 @@ int main(const int argc, const char* argv[]) RosenbrockSolver solver{ chemical_system, reactions, solver_parameters }; - auto jit{ micm::JitCompiler::Create() }; - auto start = std::chrono::high_resolution_clock::now(); JitRosenbrockSolver< GroupVectorMatrix, GroupSparseVectorMatrix, JitLinearSolver, JitProcessSet> - jit_solver(jit.get(), chemical_system, reactions, solver_parameters); + jit_solver(chemical_system, reactions, solver_parameters); auto end = std::chrono::high_resolution_clock::now(); auto jit_compile_time = std::chrono::duration_cast(end - start); diff --git a/test/tutorial/test_vectorized_matrix_solver.cpp b/test/tutorial/test_vectorized_matrix_solver.cpp index 32ba411fc..06ffade4e 100644 --- a/test/tutorial/test_vectorized_matrix_solver.cpp +++ b/test/tutorial/test_vectorized_matrix_solver.cpp @@ -15,8 +15,7 @@ using namespace micm; template using Group3Matrix = micm::VectorMatrix; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; void solve(auto& solver, auto& state) { diff --git a/test/unit/jit/test_jit_function.cpp b/test/unit/jit/test_jit_function.cpp index b1fe120b9..9e695af24 100644 --- a/test/unit/jit/test_jit_function.cpp +++ b/test/unit/jit/test_jit_function.cpp @@ -6,13 +6,7 @@ // This test creates a function that adds two integers and returns the sum TEST(JitFunction, SimpleInt32Function) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_int32") .SetArguments({ { "foo", micm::JitType::Int32 }, { "bar", micm::JitType::Int32 } }) .SetReturnType(micm::JitType::Int32); @@ -30,13 +24,7 @@ TEST(JitFunction, SimpleInt32Function) // This test creates a function that adds two integers and returns the sum TEST(JitFunction, SimpleInt64Function) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_int64") .SetArguments({ { "foo", micm::JitType::Int64 }, { "bar", micm::JitType::Int64 } }) .SetReturnType(micm::JitType::Int64); @@ -53,13 +41,7 @@ TEST(JitFunction, SimpleInt64Function) // This test creates a function that adds two floats and returns the sum TEST(JitFunction, SimpleFloatFunction) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_float") .SetArguments({ { "foo", micm::JitType::Float }, { "bar", micm::JitType::Float } }) .SetReturnType(micm::JitType::Float); @@ -76,13 +58,7 @@ TEST(JitFunction, SimpleFloatFunction) // This test creates a function that adds two doubles and returns the sum TEST(JitFunction, SimpleDoubleFunction) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_double") .SetArguments({ { "foo", micm::JitType::Double }, { "bar", micm::JitType::Double } }) .SetReturnType(micm::JitType::Double); @@ -99,13 +75,7 @@ TEST(JitFunction, SimpleDoubleFunction) // This test creates a function that returns an OR of two booleans TEST(JitFunction, SimpleBoolFunction) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_bool") .SetArguments({ { "foo", micm::JitType::Bool }, { "bar", micm::JitType::Bool } }) .SetReturnType(micm::JitType::Bool); @@ -123,13 +93,7 @@ TEST(JitFunction, SimpleBoolFunction) // second element of the second array as the sum, and also returns the sum TEST(JitFunction, SimpleInt32PtrFunction) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_int32_ptr") .SetArguments({ { "foo", micm::JitType::Int32Ptr }, { "bar", micm::JitType::Int32Ptr } }) .SetReturnType(micm::JitType::Int32); @@ -160,13 +124,7 @@ TEST(JitFunction, SimpleInt32PtrFunction) // second element of the second array as the sum, and also returns the sum TEST(JitFunction, SimpleInt64PtrFunction) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_int64_ptr") .SetArguments({ { "foo", micm::JitType::Int64Ptr }, { "bar", micm::JitType::Int64Ptr } }) .SetReturnType(micm::JitType::Int64); @@ -197,13 +155,7 @@ TEST(JitFunction, SimpleInt64PtrFunction) // second element of the second array as the sum, and also returns the sum TEST(JitFunction, SimpleFloatPtrFunction) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_float_ptr") .SetArguments({ { "foo", micm::JitType::FloatPtr }, { "bar", micm::JitType::FloatPtr } }) .SetReturnType(micm::JitType::Float); @@ -234,13 +186,7 @@ TEST(JitFunction, SimpleFloatPtrFunction) // second element of the second array as the sum, and also returns the sum TEST(JitFunction, SimpleDoublePtrFunction) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo_double_ptr") .SetArguments({ { "foo", micm::JitType::DoublePtr }, { "bar", micm::JitType::DoublePtr } }) .SetReturnType(micm::JitType::Double); @@ -271,14 +217,8 @@ TEST(JitFunction, SimpleDoublePtrFunction) // and iterates 10 times and returns the summed value TEST(JitFunction, SimpleLoop) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } micm::JitFunction func = - micm::JitFunction::Create(jit.get()).SetName("foo_loop").SetArguments({}).SetReturnType(micm::JitType::Int32); + micm::JitFunction::Create().SetName("foo_loop").SetArguments({}).SetReturnType(micm::JitType::Int32); auto loop = func.StartLoop("foo loop", 0, 10); llvm::PHINode *ret_val = func.builder_->CreatePHI(func.GetType(micm::JitType::Int32), 2, "ret val"); ret_val->addIncoming(llvm::ConstantInt::get(*(func.context_), llvm::APInt(32, 1)), func.entry_block_); @@ -297,13 +237,7 @@ TEST(JitFunction, SimpleLoop) // arguments and 10 and 100, respectively TEST(JitFunction, MultipleFunctions) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction foo_func = micm::JitFunction::Create(jit.get()) + micm::JitFunction foo_func = micm::JitFunction::Create() .SetName("foo") .SetArguments({ { "arg", micm::JitType::Int32 } }) .SetReturnType(micm::JitType::Int32); @@ -312,7 +246,7 @@ TEST(JitFunction, MultipleFunctions) foo_func.builder_->CreateRet(foo_ret_val); auto foo_target = foo_func.Generate(); int32_t (*foo_func_ptr)(int) = (int32_t(*)(int))(intptr_t)foo_target.second; - micm::JitFunction bar_func = micm::JitFunction::Create(jit.get()) + micm::JitFunction bar_func = micm::JitFunction::Create() .SetName("bar") .SetArguments({ { "arg", micm::JitType::Int32 } }) .SetReturnType(micm::JitType::Int32); @@ -332,13 +266,7 @@ TEST(JitFunction, MultipleFunctions) // and returns the sum TEST(JitFunction, LocalArray) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func = micm::JitFunction::Create(jit.get()) + micm::JitFunction func = micm::JitFunction::Create() .SetName("foo") .SetArguments({ { "arg", micm::JitType::Int64 } }) .SetReturnType(micm::JitType::Int64); @@ -385,13 +313,7 @@ TEST(JitFunction, LocalArray) // add a unique number to the function argument and return the sum TEST(JitFunction, SameNameFunctions) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error] "); - EXPECT_TRUE(false); - } - micm::JitFunction func1 = micm::JitFunction::Create(jit.get()) + micm::JitFunction func1 = micm::JitFunction::Create() .SetName("foobar") .SetArguments({ { "foo", micm::JitType::Int32 } }) .SetReturnType(micm::JitType::Int32); @@ -402,7 +324,7 @@ TEST(JitFunction, SameNameFunctions) auto func1_target = func1.Generate(); int32_t (*func1_ptr)(int32_t) = (int32_t(*)(int32_t))(intptr_t)func1_target.second; - micm::JitFunction func2 = micm::JitFunction::Create(jit.get()) + micm::JitFunction func2 = micm::JitFunction::Create() .SetName("foobar") .SetArguments({ { "foo", micm::JitType::Int32 } }) .SetReturnType(micm::JitType::Int32); @@ -413,7 +335,7 @@ TEST(JitFunction, SameNameFunctions) auto func2_target = func2.Generate(); int32_t (*func2_ptr)(int32_t) = (int32_t(*)(int32_t))(intptr_t)func2_target.second; - micm::JitFunction func3 = micm::JitFunction::Create(jit.get()) + micm::JitFunction func3 = micm::JitFunction::Create() .SetName("foobar") .SetArguments({ { "foo", micm::JitType::Int32 } }) .SetReturnType(micm::JitType::Int32); diff --git a/test/unit/openmp/test_openmp_jit_solver.cpp b/test/unit/openmp/test_openmp_jit_solver.cpp index 2fbee6cc8..2cc7acc3e 100644 --- a/test/unit/openmp/test_openmp_jit_solver.cpp +++ b/test/unit/openmp/test_openmp_jit_solver.cpp @@ -12,8 +12,7 @@ using namespace micm; template using Group1VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; TEST(OpenMP, JITOneSolverManyStates) { @@ -22,11 +21,7 @@ TEST(OpenMP, JITOneSolverManyStates) SolverConfig solverConfig; std::string config_path = "./unit_configs/robertson"; - ConfigParseStatus status = solverConfig.ReadAndParse(config_path); - if (status != micm::ConfigParseStatus::Success) - { - throw "Parsing failed"; - } + solverConfig.ReadAndParse(config_path); micm::SolverParameters solver_params = solverConfig.GetSolverParams(); @@ -35,13 +30,12 @@ TEST(OpenMP, JITOneSolverManyStates) std::vector> results(n_threads); - auto jit{ micm::JitCompiler::Create() }; JitRosenbrockSolver< Group1VectorMatrix, Group1SparseVectorMatrix, JitLinearSolver<1, Group1SparseVectorMatrix>, micm::JitProcessSet<1>> - solver(jit.get(), chemical_system, reactions, RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver(chemical_system, reactions, RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); #pragma omp parallel num_threads(n_threads) { diff --git a/test/unit/process/test_cuda_process_set.cpp b/test/unit/process/test_cuda_process_set.cpp index ba7cdb8e9..ee0b84844 100644 --- a/test/unit/process/test_cuda_process_set.cpp +++ b/test/unit/process/test_cuda_process_set.cpp @@ -23,7 +23,7 @@ void compare_pair(const index_pair& a, const index_pair& b) EXPECT_EQ(a.second, b.second); } -template class CPUMatrixPolicy, template class GPUMatrixPolicy> +template void testRandomSystemAddForcingTerms(std::size_t n_cells, std::size_t n_reactions, std::size_t n_species) { auto get_n_react = std::bind(std::uniform_int_distribution<>(0, 3), std::default_random_engine()); @@ -79,15 +79,15 @@ void testRandomSystemAddForcingTerms(std::size_t n_cells, std::size_t n_reaction gpu_state.variables_.AsVector().assign(cpu_state_vars.begin(), cpu_state_vars.end()); gpu_state.variables_.CopyToDevice(); - CPUMatrixPolicy cpu_rate_constants{ n_cells, n_reactions }; - GPUMatrixPolicy gpu_rate_constants{ n_cells, n_reactions }; + CPUMatrixPolicy cpu_rate_constants{ n_cells, n_reactions }; + GPUMatrixPolicy gpu_rate_constants{ n_cells, n_reactions }; auto& cpu_rate_vars = cpu_rate_constants.AsVector(); std::ranges::generate(cpu_rate_vars, get_double); gpu_rate_constants.AsVector().assign(cpu_rate_vars.begin(), cpu_rate_vars.end()); gpu_rate_constants.CopyToDevice(); - CPUMatrixPolicy cpu_forcing{ n_cells, n_species, 1000.0 }; - GPUMatrixPolicy gpu_forcing{ n_cells, n_species, 1000.0 }; + CPUMatrixPolicy cpu_forcing{ n_cells, n_species, 1000.0 }; + GPUMatrixPolicy gpu_forcing{ n_cells, n_species, 1000.0 }; gpu_forcing.CopyToDevice(); // kernel function call @@ -110,13 +110,9 @@ void testRandomSystemAddForcingTerms(std::size_t n_cells, std::size_t n_reaction } template< - template class CPUMatrixPolicy, - template class CPUSparseMatrixPolicy, - template class GPUDenseMatrixPolicy, - template class GPUSparseMatrixPolicy> void testRandomSystemSubtractJacobianTerms(std::size_t n_cells, std::size_t n_reactions, std::size_t n_species) { @@ -174,8 +170,8 @@ void testRandomSystemSubtractJacobianTerms(std::size_t n_cells, std::size_t n_re gpu_state.variables_.AsVector().assign(cpu_state_vars.begin(), cpu_state_vars.end()); gpu_state.variables_.CopyToDevice(); - CPUMatrixPolicy cpu_rate_constants{ n_cells, n_reactions }; - GPUDenseMatrixPolicy gpu_rate_constants{ n_cells, n_reactions }; + CPUMatrixPolicy cpu_rate_constants{ n_cells, n_reactions }; + GPUDenseMatrixPolicy gpu_rate_constants{ n_cells, n_reactions }; auto& cpu_rate_vars = cpu_rate_constants.AsVector(); std::ranges::generate(cpu_rate_vars, get_double); gpu_rate_constants.AsVector().assign(cpu_rate_vars.begin(), cpu_rate_vars.end()); @@ -183,14 +179,14 @@ void testRandomSystemSubtractJacobianTerms(std::size_t n_cells, std::size_t n_re auto non_zero_elements = cpu_set.NonZeroJacobianElements(); - auto cpu_builder = CPUSparseMatrixPolicy::Create(n_species).SetNumberOfBlocks(n_cells).InitialValue(100.0); + auto cpu_builder = CPUSparseMatrixPolicy::Create(n_species).SetNumberOfBlocks(n_cells).InitialValue(100.0); for (auto& elem : non_zero_elements) cpu_builder = cpu_builder.WithElement(elem.first, elem.second); - CPUSparseMatrixPolicy cpu_jacobian{ cpu_builder }; - auto gpu_builder = GPUSparseMatrixPolicy::Create(n_species).SetNumberOfBlocks(n_cells).InitialValue(100.0); + CPUSparseMatrixPolicy cpu_jacobian{ cpu_builder }; + auto gpu_builder = GPUSparseMatrixPolicy::Create(n_species).SetNumberOfBlocks(n_cells).InitialValue(100.0); for (auto& elem : non_zero_elements) gpu_builder = gpu_builder.WithElement(elem.first, elem.second); - GPUSparseMatrixPolicy gpu_jacobian{ gpu_builder }; + GPUSparseMatrixPolicy gpu_jacobian{ gpu_builder }; gpu_jacobian.CopyToDevice(); cpu_set.SetJacobianFlatIds(cpu_jacobian); @@ -212,17 +208,10 @@ void testRandomSystemSubtractJacobianTerms(std::size_t n_cells, std::size_t n_re } } -template -using Group10000VectorMatrix = micm::VectorMatrix; - -template -using Group10000SparseVectorMatrix = micm::SparseMatrix>; - -template -using Group10000CudaDenseMatrix = micm::CudaDenseMatrix; - -template -using Group10000CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group10000VectorMatrix = micm::VectorMatrix; +using Group10000SparseVectorMatrix = micm::SparseMatrix>; +using Group10000CudaDenseMatrix = micm::CudaDenseMatrix; +using Group10000CudaSparseMatrix = micm::CudaSparseMatrix>; TEST(RandomCudaProcessSet, Forcing) { diff --git a/test/unit/process/test_jit_forcing_calculation.cpp b/test/unit/process/test_jit_forcing_calculation.cpp index be6129f92..9ebe975ce 100644 --- a/test/unit/process/test_jit_forcing_calculation.cpp +++ b/test/unit/process/test_jit_forcing_calculation.cpp @@ -17,10 +17,8 @@ #define NUM_ITERATIONS 100 -template -using ForcingTestVectorMatrix = micm::VectorMatrix; -template -using ForcingTestSparseVectorMatrix = micm::SparseMatrix>; +using ForcingTestVectorMatrix = micm::VectorMatrix; +using ForcingTestSparseVectorMatrix = micm::SparseMatrix>; TEST(JitProcessSet, ForcingFunction) { @@ -51,17 +49,11 @@ TEST(JitProcessSet, ForcingFunction) std::vector processes{ r1, r2, r3 }; - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - micm::JitProcessSet process_set{ jit.get(), processes, state.variable_map_ }; + micm::JitProcessSet process_set{processes, state.variable_map_ }; - ForcingTestVectorMatrix rate_constants{ NUM_GRID_CELLS, 3 }; - ForcingTestVectorMatrix forcing{ NUM_GRID_CELLS, 6, 0.0 }; - ForcingTestVectorMatrix jit_forcing{ NUM_GRID_CELLS, 6, 0.0 }; + ForcingTestVectorMatrix rate_constants{ NUM_GRID_CELLS, 3 }; + ForcingTestVectorMatrix forcing{ NUM_GRID_CELLS, 6, 0.0 }; + ForcingTestVectorMatrix jit_forcing{ NUM_GRID_CELLS, 6, 0.0 }; for (int i = 0; i < 3 * NUM_GRID_CELLS; ++i) rate_constants.AsVector()[i] = i * 2.0; for (int i = 0; i < 6 * NUM_GRID_CELLS; ++i) diff --git a/test/unit/process/test_jit_process_set.cpp b/test/unit/process/test_jit_process_set.cpp index 3f2e43864..37727d417 100644 --- a/test/unit/process/test_jit_process_set.cpp +++ b/test/unit/process/test_jit_process_set.cpp @@ -9,71 +9,22 @@ #include -template -using Group2VectorMatrix = micm::VectorMatrix; -template -using Group200VectorMatrix = micm::VectorMatrix; -template -using Group300VectorMatrix = micm::VectorMatrix; -template -using Group400VectorMatrix = micm::VectorMatrix; +using Group2VectorMatrix = micm::VectorMatrix; +using Group200VectorMatrix = micm::VectorMatrix; +using Group300VectorMatrix = micm::VectorMatrix; +using Group400VectorMatrix = micm::VectorMatrix; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; TEST(JitProcessSet, VectorMatrix) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testProcessSet>( - [&](const std::vector& processes, - const micm::State& state) -> micm::JitProcessSet<2> { - return micm::JitProcessSet<2>{ jit.get(), processes, state.variable_map_ }; - }); + testProcessSet>(); } TEST(RandomJitProcessSet, VectorMatrix) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testRandomSystem>( - 200, - 20, - 30, - [&](const std::vector& processes, - const micm::State& state) -> micm::JitProcessSet<200> { - return micm::JitProcessSet<200>{ jit.get(), processes, state.variable_map_ }; - }); - testRandomSystem>( - 300, - 50, - 40, - [&](const std::vector& processes, - const micm::State& state) -> micm::JitProcessSet<300> { - return micm::JitProcessSet<300>{ jit.get(), processes, state.variable_map_ }; - }); - testRandomSystem>( - 300, - 30, - 20, - [&](const std::vector& processes, - const micm::State& state) -> micm::JitProcessSet<300> { - return micm::JitProcessSet<300>{ jit.get(), processes, state.variable_map_ }; - }); - testRandomSystem>( - 400, - 100, - 80, - [&](const std::vector& processes, - const micm::State& state) -> micm::JitProcessSet<400> { - return micm::JitProcessSet<400>{ jit.get(), processes, state.variable_map_ }; - }); + testRandomSystem>(200, 20, 30); + testRandomSystem>(300, 50, 40); + testRandomSystem>(300, 30, 20); + testRandomSystem>(400, 100, 80); } \ No newline at end of file diff --git a/test/unit/process/test_process.cpp b/test/unit/process/test_process.cpp index 1061846a6..b4a7e4f01 100644 --- a/test/unit/process/test_process.cpp +++ b/test/unit/process/test_process.cpp @@ -9,7 +9,7 @@ #include -template class MatrixPolicy> +template void testProcessUpdateState(const std::size_t number_of_grid_cells) { micm::Species foo("foo", { { "molecular weight [kg mol-1]", 0.025 }, { "diffusion coefficient [m2 s-1]", 2.3e2 } }); @@ -29,14 +29,14 @@ void testProcessUpdateState(const std::size_t number_of_grid_cells) for (const auto& process : processes) for (auto& label : process.rate_constant_->CustomParameters()) param_labels.push_back(label); - micm::State state{ micm::StateParameters{ + micm::State state{ micm::StateParameters{ .number_of_grid_cells_ = number_of_grid_cells, .number_of_rate_constants_ = processes.size(), .variable_names_ = { "foo", "bar'" }, .custom_rate_parameter_labels_ = param_labels, } }; - MatrixPolicy expected_rate_constants(number_of_grid_cells, 3, 0.0); + DenseMatrixPolicy expected_rate_constants(number_of_grid_cells, 3, 0.0); std::vector params = { 0.0, 0.0, 0.0 }; auto get_double = std::bind(std::lognormal_distribution(0.0, 0.01), std::default_random_engine()); @@ -82,15 +82,15 @@ using Group4VectorMatrix = micm::VectorMatrix; TEST(Process, Matrix) { - testProcessUpdateState(5); + testProcessUpdateState>(5); } TEST(Process, VectorMatrix) { - testProcessUpdateState(5); - testProcessUpdateState(5); - testProcessUpdateState(5); - testProcessUpdateState(5); + testProcessUpdateState>(5); + testProcessUpdateState>(5); + testProcessUpdateState>(5); + testProcessUpdateState>(5); } TEST(Process, SurfaceRateConstantOnlyHasOneReactant) diff --git a/test/unit/process/test_process_set.cpp b/test/unit/process/test_process_set.cpp index 4a9973094..86cc463e4 100644 --- a/test/unit/process/test_process_set.cpp +++ b/test/unit/process/test_process_set.cpp @@ -11,84 +11,34 @@ #include -template -using SparseMatrixTest = micm::SparseMatrix; +using SparseMatrixTest = micm::SparseMatrix; -template -using Group1VectorMatrix = micm::VectorMatrix; -template -using Group2VectorMatrix = micm::VectorMatrix; -template -using Group3VectorMatrix = micm::VectorMatrix; -template -using Group4VectorMatrix = micm::VectorMatrix; +using Group1VectorMatrix = micm::VectorMatrix; +using Group2VectorMatrix = micm::VectorMatrix; +using Group3VectorMatrix = micm::VectorMatrix; +using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; TEST(ProcessSet, Matrix) { - testProcessSet( - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); + testProcessSet, SparseMatrixTest, micm::ProcessSet>(); } TEST(ProcessSet, VectorMatrix) { - testProcessSet( - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); - testProcessSet( - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); - testProcessSet( - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); - testProcessSet( - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); + testProcessSet(); + testProcessSet(); + testProcessSet(); + testProcessSet(); } TEST(RandomProcessSet, Matrix) { - testRandomSystem( - 200, - 50, - 40, - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); - testRandomSystem( - 300, - 30, - 20, - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); - testRandomSystem( - 400, - 100, - 80, - [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { - return micm::ProcessSet{ processes, state.variable_map_ }; - }); + testRandomSystem, SparseMatrixTest, micm::ProcessSet>(200, 50, 40); + testRandomSystem, SparseMatrixTest, micm::ProcessSet>(300, 30, 20); + testRandomSystem, SparseMatrixTest, micm::ProcessSet>(400, 100, 80); } \ No newline at end of file diff --git a/test/unit/process/test_process_set_policy.hpp b/test/unit/process/test_process_set_policy.hpp index 155b18f33..1139686f3 100644 --- a/test/unit/process/test_process_set_policy.hpp +++ b/test/unit/process/test_process_set_policy.hpp @@ -13,10 +13,8 @@ void compare_pair(const index_pair& a, const index_pair& b) EXPECT_EQ(a.second, b.second); } -template class MatrixPolicy, template class SparseMatrixPolicy, class ProcessSetPolicy> -void testProcessSet(const std::function&, - const micm::State&)> create_set) +template +void testProcessSet() { auto foo = micm::Species("foo"); auto bar = micm::Species("bar"); @@ -29,7 +27,7 @@ void testProcessSet(const std::function{ foo, bar, qux, baz, quz, quuz, corge } }; - micm::State state( + micm::State state( micm::StateParameters{ .number_of_grid_cells_ = 2, .number_of_rate_constants_ = 3, .variable_names_{ "foo", "bar", "baz", "quz", "quuz", "corge" } }); @@ -57,19 +55,19 @@ void testProcessSet(const std::function{ r1, r2, r3 }, state); + ProcessSetPolicy set = ProcessSetPolicy(std::vector{ r1, r2, r3 }, state.variable_map_); EXPECT_EQ(state.variables_.NumRows(), 2); EXPECT_EQ(state.variables_.NumColumns(), 6); state.variables_[0] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.0 }; state.variables_[1] = { 1.1, 1.2, 1.3, 1.4, 1.5, 0.0 }; - MatrixPolicy rate_constants{ 2, 3 }; + DenseMatrixPolicy rate_constants{ 2, 3 }; rate_constants[0] = { 10.0, 20.0, 30.0 }; rate_constants[1] = { 110.0, 120.0, 130.0 }; - MatrixPolicy forcing{ 2, 5, 1000.0 }; + DenseMatrixPolicy forcing{ 2, 5, 1000.0 }; - set.template AddForcingTerms(rate_constants, state.variables_, forcing); + set.template AddForcingTerms(rate_constants, state.variables_, forcing); EXPECT_EQ(forcing[0][0], 1000.0 - 10.0 * 0.1 * 0.3 + 20.0 * 0.2); EXPECT_EQ(forcing[1][0], 1000.0 - 110.0 * 1.1 * 1.3 + 120.0 * 1.2); EXPECT_EQ(forcing[0][1], 1000.0 + 10.0 * 0.1 * 0.3 - 20.0 * 0.2); @@ -103,10 +101,10 @@ void testProcessSet(const std::function::Create(5).SetNumberOfBlocks(2).InitialValue(100.0); + auto builder = SparseMatrixPolicy::Create(5).SetNumberOfBlocks(2).InitialValue(100.0); for (auto& elem : non_zero_elements) builder = builder.WithElement(elem.first, elem.second); - SparseMatrixPolicy jacobian{ builder }; + SparseMatrixPolicy jacobian{ builder }; set.SetJacobianFlatIds(jacobian); set.SubtractJacobianTerms(rate_constants, state.variables_, jacobian); EXPECT_DOUBLE_EQ(jacobian[0][0][0], 100.0 + 10.0 * 0.3); // foo -> foo @@ -135,14 +133,8 @@ void testProcessSet(const std::function class MatrixPolicy, template class SparseMatrixPolicy, class ProcessSetPolicy> -void testRandomSystem( - std::size_t n_cells, - std::size_t n_reactions, - std::size_t n_species, - const std::function< - ProcessSetPolicy(const std::vector&, const micm::State&)> - create_set) +template +void testRandomSystem(std::size_t n_cells, std::size_t n_reactions, std::size_t n_species) { auto get_n_react = std::bind(std::uniform_int_distribution<>(0, 3), std::default_random_engine()); auto get_n_product = std::bind(std::uniform_int_distribution<>(0, 10), std::default_random_engine()); @@ -157,7 +149,7 @@ void testRandomSystem( species_names.push_back(std::to_string(i)); } micm::Phase gas_phase{ species }; - micm::State state{ micm::StateParameters{ + micm::State state{ micm::StateParameters{ .number_of_grid_cells_ = n_cells, .number_of_rate_constants_ = n_reactions, .variable_names_{ species_names }, @@ -180,15 +172,15 @@ void testRandomSystem( auto proc = micm::Process(micm::Process::Create().SetReactants(reactants).SetProducts(products).SetPhase(gas_phase)); processes.push_back(proc); } - ProcessSetPolicy set = create_set(processes, state); + ProcessSetPolicy set = ProcessSetPolicy(processes, state.variable_map_); for (auto& elem : state.variables_.AsVector()) elem = get_double(); - MatrixPolicy rate_constants{ n_cells, n_reactions }; + DenseMatrixPolicy rate_constants{ n_cells, n_reactions }; for (auto& elem : rate_constants.AsVector()) elem = get_double(); - MatrixPolicy forcing{ n_cells, n_species, 1000.0 }; + DenseMatrixPolicy forcing{ n_cells, n_species, 1000.0 }; - set.template AddForcingTerms(rate_constants, state.variables_, forcing); + set.template AddForcingTerms(rate_constants, state.variables_, forcing); } diff --git a/test/unit/solver/CMakeLists.txt b/test/unit/solver/CMakeLists.txt index f2cc988b3..254a24b17 100644 --- a/test/unit/solver/CMakeLists.txt +++ b/test/unit/solver/CMakeLists.txt @@ -12,6 +12,7 @@ create_standard_test(NAME lu_decomposition SOURCES test_lu_decomposition.cpp) create_standard_test(NAME rosenbrock SOURCES test_rosenbrock.cpp) create_standard_test(NAME state SOURCES test_state.cpp) create_standard_test(NAME backward_euler SOURCES test_backward_euler.cpp) +create_standard_test(NAME solver_builder SOURCES test_solver_builder.cpp) # GPU tests if(MICM_ENABLE_CUDA) diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index bdf06bafd..e85cd9113 100644 --- a/test/unit/solver/test_backward_euler.cpp +++ b/test/unit/solver/test_backward_euler.cpp @@ -1,21 +1,11 @@ -#include -#include +#include #include #include #include -TEST(BackwardEuler, DefaultConstructor) -{ - EXPECT_NO_THROW(micm::BackwardEuler be); -} - -TEST(BackwardEuler, CanCallSolve) -{ - micm::BackwardEuler be; - double time_step = 1.0; - +namespace { auto a = micm::Species("A"); auto b = micm::Species("B"); auto c = micm::Species("C"); @@ -35,24 +25,29 @@ TEST(BackwardEuler, CanCallSolve) micm::ArrheniusRateConstantParameters{ .A_ = 3.3e-11, .B_ = 0, .C_ = 55 })) .SetPhase(gas_phase); - auto options = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - options.ignore_unused_species_ = true; + auto the_system = micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }); + std::vector reactions = { r1, r2 }; +} - micm::RosenbrockSolver<> rosenbrock{ micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), - std::vector{ r1, r2 }, - options }; +TEST(BackwardEuler, CanCallSolve) +{ + auto params = micm::BackwardEulerSolverParameters(); + params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; + + auto be = micm::CpuSolverBuilder, micm::SparseMatrix>() + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(1) + .SetSolverParameters(params) + .Build(); + double time_step = 1.0; - auto state = rosenbrock.GetState(); + auto state = be.GetState(); state.variables_[0] = { 1.0, 0.0, 0.0 }; state.conditions_[0].temperature_ = 272.5; state.conditions_[0].pressure_ = 101253.3; state.conditions_[0].air_density_ = 1e6; - auto linear_solver = rosenbrock.linear_solver_; - auto process_set = rosenbrock.process_set_; - auto processes = std::vector{ r1, r2 }; - - EXPECT_NO_THROW(be.Solve( - time_step, state, linear_solver, process_set, processes, rosenbrock.state_parameters_.jacobian_diagonal_elements_)); + EXPECT_NO_THROW(be.Solve(time_step, state)); } diff --git a/test/unit/solver/test_chapman_ode_solver.cpp b/test/unit/solver/test_chapman_ode_solver.cpp index 07740598e..c27cc0437 100644 --- a/test/unit/solver/test_chapman_ode_solver.cpp +++ b/test/unit/solver/test_chapman_ode_solver.cpp @@ -210,7 +210,7 @@ TEST(ChapmanMechanismHardCodedAndGeneral, dforce_dy_times_vector) void TestSolve(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.variables_[0] = { 1, 3.92e-1, 1.69e-2, 0, 3.29e1, 0, 0, 8.84, 0 }; //"M" "Ar" "CO2", "H2O", "N2", "O1D", "O", "O2", "O3", state.conditions_[0].temperature_ = 273.15; @@ -245,7 +245,7 @@ TEST(ChapmanMechanismHardCodedAndGeneral, Solve) void TestSolve10TimesLarger(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.variables_[0] = { 1, 3.92e-1, 1.69e-2, 0, 3.29e1, 0, 0, 8.84, 0 }; //"M" "Ar" "CO2", "H2O", "N2", "O1D", "O", "O2", "O3", state.conditions_[0].temperature_ = 273.15; @@ -283,7 +283,7 @@ TEST(ChapmanMechanismHardCodedAndGeneral, solve_10_times_larger) void TestSolve10TimesSmaller(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.variables_[0] = { 1, 3.92e-1, 1.69e-2, 0, 3.29e1, 0, 0, 8.84, 0 }; //"M" "Ar" "CO2", "H2O", "N2", "O1D", "O", "O2", "O3", state.conditions_[0].temperature_ = 273.15; @@ -321,7 +321,7 @@ TEST(RegressionChapmanODESolver, solve_10_times_smaller) void TestSolveWithRandomNumberDensities(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.conditions_[0].temperature_ = 273.15; state.conditions_[0].pressure_ = 1000 * 100; // 1000 hPa state.conditions_[0].air_density_ = 2.7e19; diff --git a/test/unit/solver/test_cuda_linear_solver.cpp b/test/unit/solver/test_cuda_linear_solver.cpp index 47cc56a80..eca79c22d 100644 --- a/test/unit/solver/test_cuda_linear_solver.cpp +++ b/test/unit/solver/test_cuda_linear_solver.cpp @@ -1,4 +1,3 @@ -#pragma once #include "test_linear_solver_policy.hpp" #include @@ -17,51 +16,39 @@ #include #include -template -using Group10000VectorMatrix = micm::VectorMatrix; - -template -using Group1CudaDenseMatrix = micm::CudaDenseMatrix; -template -using Group20CudaDenseMatrix = micm::CudaDenseMatrix; -template -using Group300CudaDenseMatrix = micm::CudaDenseMatrix; -template -using Group4000CudaDenseMatrix = micm::CudaDenseMatrix; -template -using Group10000CudaDenseMatrix = micm::CudaDenseMatrix; - -template -using Group10000SparseVectorMatrix = micm::SparseMatrix>; - -template -using Group1CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group20CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group300CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group4000CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group10000CudaSparseMatrix = micm::CudaSparseMatrix>; - -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> -std::vector linearSolverGenerator( - const std::function, double)> create_linear_solver, - std::size_t number_of_blocks) +using FloatingPointType = double; + +using Group10000VectorMatrix = micm::VectorMatrix; + +using Group1CudaDenseMatrix = micm::CudaDenseMatrix; +using Group20CudaDenseMatrix = micm::CudaDenseMatrix; +using Group300CudaDenseMatrix = micm::CudaDenseMatrix; +using Group4000CudaDenseMatrix = micm::CudaDenseMatrix; +using Group10000CudaDenseMatrix = micm::CudaDenseMatrix; + +using Group10000SparseVectorMatrix = micm::SparseMatrix>; + +using Group1CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group20CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group300CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group4000CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group10000CudaSparseMatrix = micm::CudaSparseMatrix>; + +template +std::vector linearSolverGenerator(std::size_t number_of_blocks) { auto gen_bool = std::bind(std::uniform_int_distribution<>(0, 1), std::default_random_engine()); auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); - auto builder = SparseMatrixPolicy::Create(10).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); + auto builder = SparseMatrixPolicy::Create(10).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) if (i == j || gen_bool()) builder = builder.WithElement(i, j); - SparseMatrixPolicy A(builder); - MatrixPolicy b(number_of_blocks, 10, 0.0); - MatrixPolicy x(number_of_blocks, 10, 100.0); + SparseMatrixPolicy A(builder); + MatrixPolicy b(number_of_blocks, 10, 0.0); + MatrixPolicy x(number_of_blocks, 10, 100.0); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) @@ -74,25 +61,25 @@ std::vector linearSolverGenerator( b[i_block][i] = get_double(); // Only copy the data to the device when it is a CudaMatrix - CopyToDeviceSparse(A); - CopyToDeviceDense(b); - CopyToDeviceDense(x); + CopyToDeviceSparse(A); + CopyToDeviceDense(b); + CopyToDeviceDense(x); - LinearSolverPolicy solver = create_linear_solver(A, 1.0e-30); - std::pair, SparseMatrixPolicy> lu = - micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); - SparseMatrixPolicy lower_matrix = std::move(lu.first); - SparseMatrixPolicy upper_matrix = std::move(lu.second); + LinearSolverPolicy solver = LinearSolverPolicy(A, 1.0e-30); + std::pair lu = + micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); + SparseMatrixPolicy lower_matrix = std::move(lu.first); + SparseMatrixPolicy upper_matrix = std::move(lu.second); // Only copy the data to the device when it is a CudaMatrix - CopyToDeviceSparse(lower_matrix); - CopyToDeviceSparse(upper_matrix); + CopyToDeviceSparse(lower_matrix); + CopyToDeviceSparse(upper_matrix); solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(b, x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix - CopyToHostDense(x); + CopyToHostDense(x); return x.AsVector(); } @@ -103,22 +90,12 @@ void verify_gpu_against_cpu() std::vector cpu_x = linearSolverGenerator< Group10000VectorMatrix, Group10000SparseVectorMatrix, - micm::LinearSolver>( - [](const Group10000SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 10000); + micm::LinearSolver>(10000); std::vector gpu_x = linearSolverGenerator< Group10000CudaDenseMatrix, Group10000CudaSparseMatrix, - micm::CudaLinearSolver>( - [](const Group10000CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 10000); + micm::CudaLinearSolver>(10000); for (int i = 0; i < cpu_x.size(); i++) { @@ -131,82 +108,23 @@ TEST(CudaLinearSolver, DenseMatrixVectorOrderingPolicy) testDenseMatrix< Group1CudaDenseMatrix, Group1CudaSparseMatrix, - micm::CudaLinearSolver>( - [](const Group1CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }); + micm::CudaLinearSolver>(); } TEST(CudaLinearSolver, RandomMatrixVectorOrderingPolicy) { - testRandomMatrix>( - [](const Group1CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 1); - testRandomMatrix>( - [](const Group20CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 20); - testRandomMatrix< - Group300CudaDenseMatrix, - Group300CudaSparseMatrix, - micm::CudaLinearSolver>( - [](const Group300CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 300); - testRandomMatrix< - Group4000CudaDenseMatrix, - Group4000CudaSparseMatrix, - micm::CudaLinearSolver>( - [](const Group4000CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 4000); + testRandomMatrix>(1); + testRandomMatrix>(20); + testRandomMatrix>(300); + testRandomMatrix>(4000); } TEST(CudaLinearSolver, DiagonalMatrixVectorOrderingPolicy) { - testDiagonalMatrix>( - [](const Group1CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 1); - testDiagonalMatrix< - Group20CudaDenseMatrix, - Group20CudaSparseMatrix, - micm::CudaLinearSolver>( - [](const Group20CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 20); - testDiagonalMatrix< - Group300CudaDenseMatrix, - Group300CudaSparseMatrix, - micm::CudaLinearSolver>( - [](const Group300CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 300); - testDiagonalMatrix< - Group4000CudaDenseMatrix, - Group4000CudaSparseMatrix, - micm::CudaLinearSolver>( - [](const Group4000CudaSparseMatrix& matrix, - double initial_value) -> micm::CudaLinearSolver { - return micm::CudaLinearSolver{ matrix, initial_value }; - }, - 4000); + testDiagonalMatrix>(1); + testDiagonalMatrix>(20); + testDiagonalMatrix>(300); + testDiagonalMatrix>(4000); } TEST(CudaLinearSolver, RandomMatrixVectorOrderingForGPU) diff --git a/test/unit/solver/test_cuda_lu_decomposition.cpp b/test/unit/solver/test_cuda_lu_decomposition.cpp index 1ac512d4d..cc7a79146 100644 --- a/test/unit/solver/test_cuda_lu_decomposition.cpp +++ b/test/unit/solver/test_cuda_lu_decomposition.cpp @@ -1,4 +1,4 @@ -#pragma once +#include "test_lu_decomposition_policy.hpp" #include #include @@ -13,59 +13,20 @@ #include #include -template class SparseMatrixPolicy> -void check_results( - const SparseMatrixPolicy& A, - const SparseMatrixPolicy& L, - const SparseMatrixPolicy& U, - const std::function f) -{ - EXPECT_EQ(A.NumberOfBlocks(), L.NumberOfBlocks()); - EXPECT_EQ(A.NumberOfBlocks(), U.NumberOfBlocks()); - for (std::size_t i_block = 0; i_block < A.NumberOfBlocks(); ++i_block) - { - for (std::size_t i = 0; i < A.NumRows(); ++i) - { - for (std::size_t j = 0; j < A.NumColumns(); ++j) - { - T result{}; - for (std::size_t k = 0; k < A.NumRows(); ++k) - { - if (!(L.IsZero(i, k) || U.IsZero(k, j))) - { - result += L[i_block][i][k] * U[i_block][k][j]; - } - } - // Make sure these are actually triangular matrices - EXPECT_TRUE(i >= j || L.IsZero(i, j)); - EXPECT_TRUE(j >= i || U.IsZero(i, j)); - if (A.IsZero(i, j)) - { - f(result, T{}); - } - else - { - f(result, A[i_block][i][j]); - } - } - } - } -} - -template class CPUSparseMatrixPolicy, template class GPUSparseMatrixPolicy> -void testRandomMatrix(size_t n_grids) +template +void testCudaRandomMatrix(size_t n_grids) { auto gen_bool = std::bind(std::uniform_int_distribution<>(0, 1), std::default_random_engine()); auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); - auto builder = CPUSparseMatrixPolicy::Create(10).SetNumberOfBlocks(n_grids).InitialValue(1.0e-30); + auto builder = CPUSparseMatrixPolicy::Create(10).SetNumberOfBlocks(n_grids).InitialValue(1.0e-30); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) if (i == j || gen_bool()) builder = builder.WithElement(i, j); - CPUSparseMatrixPolicy cpu_A(builder); - GPUSparseMatrixPolicy gpu_A(builder); + CPUSparseMatrixPolicy cpu_A(builder); + GPUSparseMatrixPolicy gpu_A(builder); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) @@ -81,15 +42,15 @@ void testRandomMatrix(size_t n_grids) gpu_A.CopyToDevice(); gpu_LU.first.CopyToDevice(); gpu_LU.second.CopyToDevice(); - gpu_lud.Decompose(gpu_A, gpu_LU.first, gpu_LU.second); + gpu_lud.Decompose(gpu_A, gpu_LU.first, gpu_LU.second); gpu_LU.first.CopyToHost(); gpu_LU.second.CopyToHost(); - check_results( + check_results( gpu_A, gpu_LU.first, gpu_LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); - micm::LuDecomposition cpu_lud = micm::LuDecomposition::Create(cpu_A); - auto cpu_LU = micm::LuDecomposition::GetLUMatrices(cpu_A, 1.0e-30); - cpu_lud.Decompose(cpu_A, cpu_LU.first, cpu_LU.second); + micm::LuDecomposition cpu_lud = micm::LuDecomposition::Create(cpu_A); + auto cpu_LU = micm::LuDecomposition::GetLUMatrices(cpu_A, 1.0e-30); + cpu_lud.Decompose(cpu_A, cpu_LU.first, cpu_LU.second); // checking GPU result again CPU size_t L_size = cpu_LU.first.AsVector().size(); @@ -108,28 +69,20 @@ void testRandomMatrix(size_t n_grids) }; } -template -using Group1CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group100CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group1000CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group100000CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group1CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group100CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group1000CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group100000CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group1CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group100CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group1000CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group100000CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group1CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group100CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group1000CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group100000CudaSparseMatrix = micm::CudaSparseMatrix>; TEST(CudaLuDecomposition, RandomMatrixVectorOrdering) { - testRandomMatrix(1); - testRandomMatrix(100); - testRandomMatrix(1000); - testRandomMatrix(100000); + testCudaRandomMatrix(1); + testCudaRandomMatrix(100); + testCudaRandomMatrix(1000); + testCudaRandomMatrix(100000); } \ No newline at end of file diff --git a/test/unit/solver/test_cuda_rosenbrock.cpp b/test/unit/solver/test_cuda_rosenbrock.cpp index 65dcdb303..0f3b92863 100644 --- a/test/unit/solver/test_cuda_rosenbrock.cpp +++ b/test/unit/solver/test_cuda_rosenbrock.cpp @@ -31,23 +31,15 @@ using Group300GPUVectorMatrix = micm::CudaDenseMatrix; template using Group4000GPUVectorMatrix = micm::CudaDenseMatrix; -template -using Group1CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group20CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group300CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group4000CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group1CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group20CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group300CPUSparseVectorMatrix = micm::SparseMatrix>; +using Group4000CPUSparseVectorMatrix = micm::SparseMatrix>; -template -using Group1GPUSparseVectorMatrix = micm::CudaSparseMatrix>; -template -using Group20GPUSparseVectorMatrix = micm::CudaSparseMatrix>; -template -using Group300GPUSparseVectorMatrix = micm::CudaSparseMatrix>; -template -using Group4000GPUSparseVectorMatrix = micm::CudaSparseMatrix>; +using Group1GPUSparseVectorMatrix = micm::CudaSparseMatrix>; +using Group20GPUSparseVectorMatrix = micm::CudaSparseMatrix>; +using Group300GPUSparseVectorMatrix = micm::CudaSparseMatrix>; +using Group4000GPUSparseVectorMatrix = micm::CudaSparseMatrix>; // the following alias works for a CudaDenserMatrix with given row and any columns template @@ -76,30 +68,18 @@ template using Group3395043CudaDenseMatrix = micm::CudaDenseMatrix; // the following alias works for a CudaSparseMatrix with given rows and any columns -template -using Group1CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group2CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group4CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group7CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group12CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group16CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group20CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group5599CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group6603CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group200041CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group421875CudaSparseMatrix = micm::CudaSparseMatrix>; -template -using Group3395043CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group1CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group2CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group4CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group7CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group12CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group16CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group20CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group5599CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group6603CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group200041CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group421875CudaSparseMatrix = micm::CudaSparseMatrix>; +using Group3395043CudaSparseMatrix = micm::CudaSparseMatrix>; template RosenbrockPolicy getSolver(std::size_t number_of_grid_cells) @@ -149,11 +129,9 @@ RosenbrockPolicy getSolver(std::size_t number_of_grid_cells) template< template class CPUMatrixPolicy, - template class CPUSparseMatrixPolicy, template class GPUMatrixPolicy, - template class GPUSparseMatrixPolicy> void testAlphaMinusJacobian(std::size_t number_of_grid_cells) { @@ -223,7 +201,7 @@ void testAlphaMinusJacobian(std::size_t number_of_grid_cells) // In this test, all the elements in the same array are identical; // thus the calculated RMSE should be the same no matter what the size of the array is. -template class MatrixPolicy, template class SparseMatrixPolicy> +template class MatrixPolicy, class SparseMatrixPolicy> void testNormalizedErrorConst(const size_t number_of_grid_cells) { auto gpu_solver = getSolver>(number_of_grid_cells); @@ -264,7 +242,7 @@ void testNormalizedErrorConst(const size_t number_of_grid_cells) // In this test, the elements in the same array are different; // thus the calculated RMSE will change when the size of the array changes. -template class MatrixPolicy, template class SparseMatrixPolicy> +template class MatrixPolicy, class SparseMatrixPolicy> void testNormalizedErrorDiff(const size_t number_of_grid_cells) { auto gpu_solver = getSolver>(number_of_grid_cells); diff --git a/test/unit/solver/test_jit_linear_solver.cpp b/test/unit/solver/test_jit_linear_solver.cpp index 73c98ea23..501152872 100644 --- a/test/unit/solver/test_jit_linear_solver.cpp +++ b/test/unit/solver/test_jit_linear_solver.cpp @@ -10,109 +10,33 @@ #include -template -using Group1VectorMatrix = micm::VectorMatrix; -template -using Group2VectorMatrix = micm::VectorMatrix; -template -using Group3VectorMatrix = micm::VectorMatrix; -template -using Group4VectorMatrix = micm::VectorMatrix; +using Group1VectorMatrix = micm::VectorMatrix; +using Group2VectorMatrix = micm::VectorMatrix; +using Group3VectorMatrix = micm::VectorMatrix; +using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; TEST(JitLinearSolver, DenseMatrixVectorOrdering) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testDenseMatrix< - Group1VectorMatrix, - Group1SparseVectorMatrix, - micm::JitLinearSolver<1, Group1SparseVectorMatrix, micm::JitLuDecomposition<1>>>( - [&](const Group1SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<1, Group1SparseVectorMatrix, micm::JitLuDecomposition<1>> - { - return micm::JitLinearSolver<1, Group1SparseVectorMatrix, micm::JitLuDecomposition<1>>{ jit.get(), - matrix, - initial_value }; - }); + testDenseMatrix>>(); } TEST(JitLinearSolver, RandomMatrixVectorOrdering) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testRandomMatrix>( - [&](const Group1SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<1, Group1SparseVectorMatrix> { - return micm::JitLinearSolver<1, Group1SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 1); - testRandomMatrix>( - [&](const Group2SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<2, Group2SparseVectorMatrix> { - return micm::JitLinearSolver<2, Group2SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 2); - testRandomMatrix>( - [&](const Group3SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<3, Group3SparseVectorMatrix> { - return micm::JitLinearSolver<3, Group3SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 3); - testRandomMatrix>( - [&](const Group4SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<4, Group4SparseVectorMatrix> { - return micm::JitLinearSolver<4, Group4SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 4); + testRandomMatrix>(1); + testRandomMatrix>(2); + testRandomMatrix>(3); + testRandomMatrix>(4); } TEST(JitLinearSolver, DiagonalMatrixVectorOrdering) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testDiagonalMatrix>( - [&](const Group1SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<1, Group1SparseVectorMatrix> { - return micm::JitLinearSolver<1, Group1SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 1); - testDiagonalMatrix>( - [&](const Group2SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<2, Group2SparseVectorMatrix> { - return micm::JitLinearSolver<2, Group2SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 2); - testDiagonalMatrix>( - [&](const Group3SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<3, Group3SparseVectorMatrix> { - return micm::JitLinearSolver<3, Group3SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 3); - testDiagonalMatrix>( - [&](const Group4SparseVectorMatrix& matrix, - double initial_value) -> micm::JitLinearSolver<4, Group4SparseVectorMatrix> { - return micm::JitLinearSolver<4, Group4SparseVectorMatrix>{ jit.get(), matrix, initial_value }; - }, - 4); + testDiagonalMatrix>(1); + testDiagonalMatrix>(2); + testDiagonalMatrix>(3); + testDiagonalMatrix>(4); } \ No newline at end of file diff --git a/test/unit/solver/test_jit_lu_decomposition.cpp b/test/unit/solver/test_jit_lu_decomposition.cpp index 594463aec..4c41cce1b 100644 --- a/test/unit/solver/test_jit_lu_decomposition.cpp +++ b/test/unit/solver/test_jit_lu_decomposition.cpp @@ -6,123 +6,39 @@ #include -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; TEST(JitLuDecomposition, DenseMatrixVectorOrdering) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testDenseMatrix>( - [&](const Group1SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<1> { - return micm::JitLuDecomposition<1>{ jit.get(), matrix }; - }); - testDenseMatrix>( - [&](const Group2SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<2> { - return micm::JitLuDecomposition<2>{ jit.get(), matrix }; - }); - testDenseMatrix>( - [&](const Group3SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<3> { - return micm::JitLuDecomposition<3>{ jit.get(), matrix }; - }); - testDenseMatrix>( - [&](const Group4SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<4> { - return micm::JitLuDecomposition<4>{ jit.get(), matrix }; - }); + testDenseMatrix>(); + testDenseMatrix>(); + testDenseMatrix>(); + testDenseMatrix>(); } TEST(JitLuDecomposition, SingularMatrixVectorOrdering) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testSingularMatrix>( - [&](const Group1SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<1> { - return micm::JitLuDecomposition<1>{ jit.get(), matrix }; - }); - testSingularMatrix>( - [&](const Group2SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<2> { - return micm::JitLuDecomposition<2>{ jit.get(), matrix }; - }); - testSingularMatrix>( - [&](const Group3SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<3> { - return micm::JitLuDecomposition<3>{ jit.get(), matrix }; - }); - testSingularMatrix>( - [&](const Group4SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<4> { - return micm::JitLuDecomposition<4>{ jit.get(), matrix }; - }); + testSingularMatrix>(); + testSingularMatrix>(); + testSingularMatrix>(); + testSingularMatrix>(); } TEST(JitLuDecomposition, RandomMatrixVectorOrdering) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testRandomMatrix>( - [&](const Group1SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<1> { - return micm::JitLuDecomposition<1>{ jit.get(), matrix }; - }, - 1); - testRandomMatrix>( - [&](const Group2SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<2> { - return micm::JitLuDecomposition<2>{ jit.get(), matrix }; - }, - 2); - testRandomMatrix>( - [&](const Group3SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<3> { - return micm::JitLuDecomposition<3>{ jit.get(), matrix }; - }, - 3); - testRandomMatrix>( - [&](const Group4SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<4> { - return micm::JitLuDecomposition<4>{ jit.get(), matrix }; - }, - 4); + testRandomMatrix>(1); + testRandomMatrix>(2); + testRandomMatrix>(3); + testRandomMatrix>(4); } TEST(JitLuDecomposition, DiagonalMatrixVectorOrdering) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testDiagonalMatrix>( - [&](const Group1SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<1> { - return micm::JitLuDecomposition<1>{ jit.get(), matrix }; - }, - 1); - testDiagonalMatrix>( - [&](const Group2SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<2> { - return micm::JitLuDecomposition<2>{ jit.get(), matrix }; - }, - 2); - testDiagonalMatrix>( - [&](const Group3SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<3> { - return micm::JitLuDecomposition<3>{ jit.get(), matrix }; - }, - 3); - testDiagonalMatrix>( - [&](const Group4SparseVectorMatrix& matrix) -> micm::JitLuDecomposition<4> { - return micm::JitLuDecomposition<4>{ jit.get(), matrix }; - }, - 4); + testDiagonalMatrix>(1); + testDiagonalMatrix>(2); + testDiagonalMatrix>(3); + testDiagonalMatrix>(4); } diff --git a/test/unit/solver/test_jit_rosenbrock.cpp b/test/unit/solver/test_jit_rosenbrock.cpp index dfb2dc6c8..b160e7887 100644 --- a/test/unit/solver/test_jit_rosenbrock.cpp +++ b/test/unit/solver/test_jit_rosenbrock.cpp @@ -10,13 +10,13 @@ #include -template class MatrixPolicy, template class SparseMatrixPolicy> +template class MatrixPolicy, class SparseMatrixPolicy> micm::JitRosenbrockSolver< MatrixPolicy, SparseMatrixPolicy, micm::JitLinearSolver, micm::JitProcessSet> -getSolver(std::shared_ptr jit) +getSolver() { // ---- foo bar baz quz quuz // foo 0 1 2 - - @@ -53,7 +53,6 @@ getSolver(std::shared_ptr jit) SparseMatrixPolicy, micm::JitLinearSolver, micm::JitProcessSet>( - jit, micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3 }, micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(number_of_grid_cells, false)); @@ -62,10 +61,10 @@ getSolver(std::shared_ptr jit) template using SparseMatrix = micm::SparseMatrix; -template class MatrixPolicy, template class SparseMatrixPolicy> -void testAlphaMinusJacobian(std::shared_ptr jit) +template class MatrixPolicy, class SparseMatrixPolicy> +void testAlphaMinusJacobian() { - auto solver = getSolver(jit); + auto solver = getSolver(); // return; auto jacobian = solver.GetState().jacobian_; @@ -121,14 +120,10 @@ using Group3VectorMatrix = micm::VectorMatrix; template using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; void run_solver(auto& solver) { @@ -157,22 +152,14 @@ void run_solver(auto& solver) TEST(JitRosenbrockSolver, AlphaMinusJacobian) { - auto jit{ micm::JitCompiler::Create() }; - if (auto err = jit.takeError()) - { - llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); - EXPECT_TRUE(false); - } - testAlphaMinusJacobian<1, Group1VectorMatrix, Group1SparseVectorMatrix>(jit.get()); - testAlphaMinusJacobian<2, Group2VectorMatrix, Group2SparseVectorMatrix>(jit.get()); - testAlphaMinusJacobian<3, Group3VectorMatrix, Group3SparseVectorMatrix>(jit.get()); - testAlphaMinusJacobian<4, Group4VectorMatrix, Group4SparseVectorMatrix>(jit.get()); + testAlphaMinusJacobian<1, Group1VectorMatrix, Group1SparseVectorMatrix>(); + testAlphaMinusJacobian<2, Group2VectorMatrix, Group2SparseVectorMatrix>(); + testAlphaMinusJacobian<3, Group3VectorMatrix, Group3SparseVectorMatrix>(); + testAlphaMinusJacobian<4, Group4VectorMatrix, Group4SparseVectorMatrix>(); } TEST(JitRosenbrockSolver, MultipleInstances) { - auto jit{ micm::JitCompiler::Create() }; - micm::SolverConfig solverConfig; std::string config_path = "./unit_configs/robertson"; solverConfig.ReadAndParse(config_path); @@ -189,19 +176,19 @@ TEST(JitRosenbrockSolver, MultipleInstances) Group1SparseVectorMatrix, micm::JitLinearSolver<1, Group1SparseVectorMatrix>, micm::JitProcessSet<1>> - solver1(jit.get(), chemical_system, reactions, solver_parameters); + solver1(chemical_system, reactions, solver_parameters); micm::JitRosenbrockSolver< Group1VectorMatrix, Group1SparseVectorMatrix, micm::JitLinearSolver<1, Group1SparseVectorMatrix>, micm::JitProcessSet<1>> - solver2(jit.get(), chemical_system, reactions, solver_parameters); + solver2(chemical_system, reactions, solver_parameters); micm::JitRosenbrockSolver< Group1VectorMatrix, Group1SparseVectorMatrix, micm::JitLinearSolver<1, Group1SparseVectorMatrix>, micm::JitProcessSet<1>> - solver3(jit.get(), chemical_system, reactions, solver_parameters); + solver3(chemical_system, reactions, solver_parameters); run_solver(solver1); run_solver(solver2); diff --git a/test/unit/solver/test_linear_solver.cpp b/test/unit/solver/test_linear_solver.cpp index 7f59bea53..c1268934c 100644 --- a/test/unit/solver/test_linear_solver.cpp +++ b/test/unit/solver/test_linear_solver.cpp @@ -10,136 +10,63 @@ #include -template -using SparseMatrixTest = micm::SparseMatrix; +using FloatingPointType = double; + +using DenseMatrixTest = micm::Matrix; +using SparseMatrixTest = micm::SparseMatrix; TEST(LinearSolver, DenseMatrixStandardOrdering) { - testDenseMatrix>( - [](const SparseMatrixTest& matrix, double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }); + testDenseMatrix>(); } TEST(LinearSolver, RandomMatrixStandardOrdering) { - testRandomMatrix>( - [](const SparseMatrixTest& matrix, double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); + testRandomMatrix>(5); } TEST(LinearSolver, DiagonalMatrixStandardOrdering) { - testDiagonalMatrix>( - [](const SparseMatrixTest& matrix, double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); + testDiagonalMatrix>(5); } TEST(LinearSolver, DiagonalMarkowitzReorder) { - testMarkowitzReordering(); + testMarkowitzReordering, SparseMatrixTest>(); } -template -using Group1VectorMatrix = micm::VectorMatrix; -template -using Group2VectorMatrix = micm::VectorMatrix; -template -using Group3VectorMatrix = micm::VectorMatrix; -template -using Group4VectorMatrix = micm::VectorMatrix; +using Group1VectorMatrix = micm::VectorMatrix; +using Group2VectorMatrix = micm::VectorMatrix; +using Group3VectorMatrix = micm::VectorMatrix; +using Group4VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; TEST(LinearSolver, DenseMatrixVectorOrdering) { - testDenseMatrix>( - [](const Group1SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }); - testDenseMatrix>( - [](const Group2SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }); - testDenseMatrix>( - [](const Group3SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }); - testDenseMatrix>( - [](const Group4SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }); + testDenseMatrix>(); + testDenseMatrix>(); + testDenseMatrix>(); + testDenseMatrix>(); } TEST(LinearSolver, RandomMatrixVectorOrdering) { - testRandomMatrix>( - [](const Group1SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); - testRandomMatrix>( - [](const Group2SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); - testRandomMatrix>( - [](const Group3SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); - testRandomMatrix>( - [](const Group4SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); + testRandomMatrix>(5); + testRandomMatrix>(5); + testRandomMatrix>(5); + testRandomMatrix>(5); } TEST(LinearSolver, DiagonalMatrixVectorOrdering) { - testDiagonalMatrix>( - [](const Group1SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); - testDiagonalMatrix>( - [](const Group2SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); - testDiagonalMatrix>( - [](const Group3SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); - testDiagonalMatrix>( - [](const Group4SparseVectorMatrix& matrix, - double initial_value) -> micm::LinearSolver { - return micm::LinearSolver{ matrix, initial_value }; - }, - 5); + testDiagonalMatrix>(5); + testDiagonalMatrix>(5); + testDiagonalMatrix>(5); + testDiagonalMatrix>(5); } TEST(LinearSolver, VectorDiagonalMarkowitzReordering) diff --git a/test/unit/solver/test_linear_solver_policy.hpp b/test/unit/solver/test_linear_solver_policy.hpp index 0adc83699..9fe8f548d 100644 --- a/test/unit/solver/test_linear_solver_policy.hpp +++ b/test/unit/solver/test_linear_solver_policy.hpp @@ -8,8 +8,8 @@ // Define the following three functions that only work for the CudaMatrix; the if constexpr statement is evalauted at // compile-time Reference: https://www.modernescpp.com/index.php/using-requires-expression-in-c-20-as-a-standalone-feature/ -template class MatrixPolicy> -void CopyToDeviceDense(MatrixPolicy& matrix) +template +void CopyToDeviceDense(MatrixPolicy& matrix) { if constexpr (requires { { @@ -19,8 +19,8 @@ void CopyToDeviceDense(MatrixPolicy& matrix) matrix.CopyToDevice(); } -template class SparseMatrixPolicy> -void CopyToDeviceSparse(SparseMatrixPolicy& matrix) +template +void CopyToDeviceSparse(SparseMatrixPolicy& matrix) { if constexpr (requires { { @@ -30,8 +30,8 @@ void CopyToDeviceSparse(SparseMatrixPolicy& matrix) matrix.CopyToDevice(); } -template class MatrixPolicy> -void CopyToHostDense(MatrixPolicy& matrix) +template +void CopyToHostDense(MatrixPolicy& matrix) { if constexpr (requires { { @@ -41,11 +41,11 @@ void CopyToHostDense(MatrixPolicy& matrix) matrix.CopyToHost(); } -template class MatrixPolicy, template class SparseMatrixPolicy> +template void check_results( - const SparseMatrixPolicy A, - const MatrixPolicy b, - const MatrixPolicy x, + const SparseMatrixPolicy A, + const MatrixPolicy b, + const MatrixPolicy x, const std::function f) { T result; @@ -64,8 +64,8 @@ void check_results( } } -template class SparseMatrixPolicy> -void print_matrix(const SparseMatrixPolicy& matrix, std::size_t width) +template +void print_matrix(const SparseMatrixPolicy& matrix, std::size_t width) { for (std::size_t i_block = 0; i_block < matrix.NumberOfBlocks(); ++i_block) { @@ -90,10 +90,12 @@ void print_matrix(const SparseMatrixPolicy& matrix, std::size_t width) } } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> -void testDenseMatrix(const std::function, double)> create_linear_solver) +template +void testDenseMatrix() { - SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(3) + using FloatingPointType = typename MatrixPolicy::value_type; + + SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(3) .InitialValue(1.0e-30) .WithElement(0, 0) .WithElement(0, 1) @@ -104,8 +106,8 @@ void testDenseMatrix(const std::function b(1, 3, 0.0); - MatrixPolicy x(1, 3, 100.0); + MatrixPolicy b(1, 3, 0.0); + MatrixPolicy x(1, 3, 100.0); A[0][0][0] = 2; A[0][0][1] = -1; @@ -122,46 +124,46 @@ void testDenseMatrix(const std::function(A); - CopyToDeviceDense(b); - CopyToDeviceDense(x); + CopyToDeviceSparse(A); + CopyToDeviceDense(b); + CopyToDeviceDense(x); - LinearSolverPolicy solver = create_linear_solver(A, 1.0e-30); - auto lu = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); + LinearSolverPolicy solver = LinearSolverPolicy(A, 1.0e-30); + auto lu = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); // Only copy the data to the device when it is a CudaMatrix - CopyToDeviceSparse(lower_matrix); - CopyToDeviceSparse(upper_matrix); + CopyToDeviceSparse(lower_matrix); + CopyToDeviceSparse(upper_matrix); solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(b, x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix - CopyToHostDense(x); + CopyToHostDense(x); - check_results( - A, b, x, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); + check_results( + A, b, x, [&](const FloatingPointType a, const FloatingPointType b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> -void testRandomMatrix( - const std::function, double)> create_linear_solver, - std::size_t number_of_blocks) +template +void testRandomMatrix(std::size_t number_of_blocks) { + using FloatingPointType = typename MatrixPolicy::value_type; + auto gen_bool = std::bind(std::uniform_int_distribution<>(0, 1), std::default_random_engine()); auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); - auto builder = SparseMatrixPolicy::Create(10).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); + auto builder = SparseMatrixPolicy::Create(10).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) if (i == j || gen_bool()) builder = builder.WithElement(i, j); - SparseMatrixPolicy A(builder); - MatrixPolicy b(number_of_blocks, 10, 0.0); - MatrixPolicy x(number_of_blocks, 10, 100.0); + SparseMatrixPolicy A(builder); + MatrixPolicy b(number_of_blocks, 10, 0.0); + MatrixPolicy x(number_of_blocks, 10, 100.0); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) @@ -174,78 +176,78 @@ void testRandomMatrix( b[i_block][i] = get_double(); // Only copy the data to the device when it is a CudaMatrix - CopyToDeviceSparse(A); - CopyToDeviceDense(b); - CopyToDeviceDense(x); + CopyToDeviceSparse(A); + CopyToDeviceDense(b); + CopyToDeviceDense(x); - LinearSolverPolicy solver = create_linear_solver(A, 1.0e-30); - auto lu = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); + LinearSolverPolicy solver = LinearSolverPolicy(A, 1.0e-30); + auto lu = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); // Only copy the data to the device when it is a CudaMatrix - CopyToDeviceSparse(lower_matrix); - CopyToDeviceSparse(upper_matrix); + CopyToDeviceSparse(lower_matrix); + CopyToDeviceSparse(upper_matrix); solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(b, x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix - CopyToHostDense(x); + CopyToHostDense(x); - check_results( - A, b, x, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); + check_results( + A, b, x, [&](const FloatingPointType a, const FloatingPointType b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); } -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> -void testDiagonalMatrix( - const std::function, double)> create_linear_solver, - std::size_t number_of_blocks) +template +void testDiagonalMatrix(std::size_t number_of_blocks) { + using FloatingPointType = typename MatrixPolicy::value_type; + auto get_double = std::bind(std::lognormal_distribution(-2.0, 4.0), std::default_random_engine()); - auto builder = SparseMatrixPolicy::Create(6).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); + auto builder = SparseMatrixPolicy::Create(6).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); for (std::size_t i = 0; i < 6; ++i) builder = builder.WithElement(i, i); - SparseMatrixPolicy A(builder); - MatrixPolicy b(number_of_blocks, 6, 0.0); - MatrixPolicy x(number_of_blocks, 6, 100.0); + SparseMatrixPolicy A(builder); + MatrixPolicy b(number_of_blocks, 6, 0.0); + MatrixPolicy x(number_of_blocks, 6, 100.0); for (std::size_t i = 0; i < 6; ++i) for (std::size_t i_block = 0; i_block < number_of_blocks; ++i_block) A[i_block][i][i] = get_double(); // Only copy the data to the device when it is a CudaMatrix - CopyToDeviceSparse(A); - CopyToDeviceDense(b); - CopyToDeviceDense(x); + CopyToDeviceSparse(A); + CopyToDeviceDense(b); + CopyToDeviceDense(x); - LinearSolverPolicy solver = create_linear_solver(A, 1.0e-30); - auto lu = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); + LinearSolverPolicy solver = LinearSolverPolicy(A, 1.0e-30); + auto lu = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); // Only copy the data to the device when it is a CudaMatrix - CopyToDeviceSparse(lower_matrix); - CopyToDeviceSparse(upper_matrix); + CopyToDeviceSparse(lower_matrix); + CopyToDeviceSparse(upper_matrix); solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(b, x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix - CopyToHostDense(x); + CopyToHostDense(x); - check_results( - A, b, x, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); + check_results( + A, b, x, [&](const FloatingPointType a, const FloatingPointType b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); } -template class MatrixPolicy, template class SparseMatrixPolicy> +template void testMarkowitzReordering() { const std::size_t order = 50; auto gen_bool = std::bind(std::uniform_int_distribution<>(0, 1), std::default_random_engine()); - MatrixPolicy orig(order, order, 0); + MatrixPolicy orig(order, order, 0); for (std::size_t i = 0; i < order; ++i) for (std::size_t j = 0; j < order; ++j) @@ -253,25 +255,25 @@ void testMarkowitzReordering() auto reorder_map = micm::DiagonalMarkowitzReorder(orig); - auto builder = SparseMatrixPolicy::Create(50); + auto builder = SparseMatrixPolicy::Create(50); for (std::size_t i = 0; i < order; ++i) for (std::size_t j = 0; j < order; ++j) if (orig[i][j] != 0) builder = builder.WithElement(i, j); - SparseMatrixPolicy orig_jac{ builder }; + SparseMatrixPolicy orig_jac{ builder }; - builder = SparseMatrixPolicy::Create(50); + builder = SparseMatrixPolicy::Create(50); for (std::size_t i = 0; i < order; ++i) for (std::size_t j = 0; j < order; ++j) if (orig[reorder_map[i]][reorder_map[j]] != 0) builder = builder.WithElement(i, j); - SparseMatrixPolicy reordered_jac{ builder }; + SparseMatrixPolicy reordered_jac{ builder }; - auto orig_LU_calc = micm::LuDecomposition::Create(orig_jac); - auto reordered_LU_calc = micm::LuDecomposition::Create(reordered_jac); + auto orig_LU_calc = micm::LuDecomposition::Create< SparseMatrixPolicy>(orig_jac); + auto reordered_LU_calc = micm::LuDecomposition::Create(reordered_jac); - auto orig_LU = orig_LU_calc.template GetLUMatrices(orig_jac, 0.0); - auto reordered_LU = reordered_LU_calc.template GetLUMatrices(reordered_jac, 0.0); + auto orig_LU = orig_LU_calc.template GetLUMatrices(orig_jac, 0.0); + auto reordered_LU = reordered_LU_calc.template GetLUMatrices(reordered_jac, 0.0); std::size_t sum_orig = 0; std::size_t sum_reordered = 0; diff --git a/test/unit/solver/test_lu_decomposition.cpp b/test/unit/solver/test_lu_decomposition.cpp index 67bdbc1b7..577ac567e 100644 --- a/test/unit/solver/test_lu_decomposition.cpp +++ b/test/unit/solver/test_lu_decomposition.cpp @@ -6,116 +6,61 @@ #include -template -using SparseMatrixTest = micm::SparseMatrix; +using SparseMatrixTest = micm::SparseMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; TEST(LuDecomposition, DenseMatrixStandardOrdering) { - testDenseMatrix( - [](const SparseMatrixTest& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); + testDenseMatrix(); } TEST(LuDecomposition, SingularMatrixStandardOrdering) { - testSingularMatrix( - [](const SparseMatrixTest& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); + testSingularMatrix(); } TEST(LuDecomposition, RandomMatrixStandardOrdering) { - testRandomMatrix( - [](const SparseMatrixTest& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); + testRandomMatrix(5); } TEST(LuDecomposition, DiagonalMatrixStandardOrdering) { - testDiagonalMatrix( - [](const SparseMatrixTest& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); + testDiagonalMatrix(5); } TEST(LuDecomposition, DenseMatrixVectorOrdering) { - testDenseMatrix( - [](const Group1SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); - testDenseMatrix( - [](const Group2SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); - testDenseMatrix( - [](const Group3SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); - testDenseMatrix( - [](const Group4SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); + testDenseMatrix(); + testDenseMatrix(); + testDenseMatrix(); + testDenseMatrix(); } TEST(LuDecomposition, SingluarMatrixVectorOrdering) { - testSingularMatrix( - [](const Group1SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); - testSingularMatrix( - [](const Group2SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); - testSingularMatrix( - [](const Group3SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); - testSingularMatrix( - [](const Group4SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }); + testSingularMatrix(); + testSingularMatrix(); + testSingularMatrix(); + testSingularMatrix(); } TEST(LuDecomposition, RandomMatrixVectorOrdering) { - testRandomMatrix( - [](const Group1SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); - testRandomMatrix( - [](const Group2SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); - testRandomMatrix( - [](const Group3SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); - testRandomMatrix( - [](const Group4SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); + testRandomMatrix(5); + testRandomMatrix(5); + testRandomMatrix(5); + testRandomMatrix(5); } TEST(LuDecomposition, DiagonalMatrixVectorOrdering) { - testDiagonalMatrix( - [](const Group1SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); - testDiagonalMatrix( - [](const Group2SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); - testDiagonalMatrix( - [](const Group3SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); - testDiagonalMatrix( - [](const Group4SparseVectorMatrix& matrix) -> micm::LuDecomposition - { return micm::LuDecomposition::Create(matrix); }, - 5); + testDiagonalMatrix(5); + testDiagonalMatrix(5); + testDiagonalMatrix(5); + testDiagonalMatrix(5); } diff --git a/test/unit/solver/test_lu_decomposition_policy.hpp b/test/unit/solver/test_lu_decomposition_policy.hpp index 9f4df81a6..eedc75c4b 100644 --- a/test/unit/solver/test_lu_decomposition_policy.hpp +++ b/test/unit/solver/test_lu_decomposition_policy.hpp @@ -7,11 +7,11 @@ #include #include -template class SparseMatrixPolicy> +template void check_results( - const SparseMatrixPolicy& A, - const SparseMatrixPolicy& L, - const SparseMatrixPolicy& U, + const SparseMatrixPolicy& A, + const SparseMatrixPolicy& L, + const SparseMatrixPolicy& U, const std::function f) { EXPECT_EQ(A.NumberOfBlocks(), L.NumberOfBlocks()); @@ -73,10 +73,10 @@ void print_matrix(const SparseMatrixPolicy& matrix, std::size_t width) } // tests example from https://www.geeksforgeeks.org/doolittle-algorithm-lu-decomposition/ -template class SparseMatrixPolicy, class LuDecompositionPolicy> -void testDenseMatrix(const std::function&)> create_lu_decomp) +template +void testDenseMatrix() { - SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(3) + SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(3) .InitialValue(1.0e-30) .WithElement(0, 0) .WithElement(0, 1) @@ -98,17 +98,17 @@ void testDenseMatrix(const std::function(A, 1.0e-30); - lud.template Decompose(A, LU.first, LU.second); + LuDecompositionPolicy lud = LuDecompositionPolicy(A); + auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); + lud.template Decompose(A, LU.first, LU.second); check_results( A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); } -template class SparseMatrixPolicy, class LuDecompositionPolicy> -void testSingularMatrix(const std::function&)> create_lu_decomp) +template +void testSingularMatrix() { - SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(2) + SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(2) .InitialValue(1.0e-30) .WithElement(0, 0) .WithElement(0, 1) @@ -120,31 +120,29 @@ void testSingularMatrix(const std::function(A, 1.0E-30); + LuDecompositionPolicy lud = LuDecompositionPolicy(A); + auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0E-30); bool is_singular{ false }; - lud.template Decompose(A, LU.first, LU.second, is_singular); + lud.template Decompose(A, LU.first, LU.second, is_singular); EXPECT_TRUE(is_singular); A[0][0][0] = 12; - lud.template Decompose(A, LU.first, LU.second, is_singular); + lud.template Decompose(A, LU.first, LU.second, is_singular); EXPECT_FALSE(is_singular); } -template class SparseMatrixPolicy, class LuDecompositionPolicy> -void testRandomMatrix( - const std::function&)> create_lu_decomp, - std::size_t number_of_blocks) +template +void testRandomMatrix(std::size_t number_of_blocks) { auto gen_bool = std::bind(std::uniform_int_distribution<>(0, 1), std::default_random_engine()); auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); - auto builder = SparseMatrixPolicy::Create(10).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); + auto builder = SparseMatrixPolicy::Create(10).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) if (i == j || gen_bool()) builder = builder.WithElement(i, j); - SparseMatrixPolicy A(builder); + SparseMatrixPolicy A(builder); for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) @@ -152,33 +150,31 @@ void testRandomMatrix( for (std::size_t i_block = 0; i_block < number_of_blocks; ++i_block) A[i_block][i][j] = get_double(); - LuDecompositionPolicy lud = create_lu_decomp(A); - auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); - lud.template Decompose(A, LU.first, LU.second); + LuDecompositionPolicy lud = LuDecompositionPolicy(A); + auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); + lud.template Decompose(A, LU.first, LU.second); check_results( A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); } -template class SparseMatrixPolicy, class LuDecompositionPolicy> -void testDiagonalMatrix( - const std::function&)> create_lu_decomp, - std::size_t number_of_blocks) +template +void testDiagonalMatrix(std::size_t number_of_blocks) { auto get_double = std::bind(std::lognormal_distribution(-2.0, 4.0), std::default_random_engine()); - auto builder = SparseMatrixPolicy::Create(6).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); + auto builder = SparseMatrixPolicy::Create(6).SetNumberOfBlocks(number_of_blocks).InitialValue(1.0e-30); for (std::size_t i = 0; i < 6; ++i) builder = builder.WithElement(i, i); - SparseMatrixPolicy A(builder); + SparseMatrixPolicy A(builder); for (std::size_t i = 0; i < 6; ++i) for (std::size_t i_block = 0; i_block < number_of_blocks; ++i_block) A[i_block][i][i] = get_double(); - LuDecompositionPolicy lud = create_lu_decomp(A); - auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); - lud.template Decompose(A, LU.first, LU.second); + LuDecompositionPolicy lud = LuDecompositionPolicy(A); + auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0e-30); + lud.template Decompose(A, LU.first, LU.second); check_results( A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); }); } \ No newline at end of file diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index a08d548ee..a4fc68a2f 100644 --- a/test/unit/solver/test_rosenbrock.cpp +++ b/test/unit/solver/test_rosenbrock.cpp @@ -7,7 +7,7 @@ #include -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> micm::RosenbrockSolver getSolver(std::size_t number_of_grid_cells) { // ---- foo bar baz quz quuz @@ -46,10 +46,9 @@ micm::RosenbrockSolver get micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(number_of_grid_cells, false)); } -template -using SparseMatrix = micm::SparseMatrix; +using SparseMatrix = micm::SparseMatrix; -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> void testAlphaMinusJacobian(std::size_t number_of_grid_cells) { auto solver = getSolver(number_of_grid_cells); @@ -105,7 +104,7 @@ void testAlphaMinusJacobian(std::size_t number_of_grid_cells) // In this test, the elements in the same array are different; // thus the calculated RMSE will change when the size of the array changes. -template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy> void testNormalizedErrorDiff(const size_t number_of_grid_cells) { auto solver = getSolver(number_of_grid_cells); @@ -149,10 +148,10 @@ void testNormalizedErrorDiff(const size_t number_of_grid_cells) TEST(RosenbrockSolver, StandardAlphaMinusJacobian) { - testAlphaMinusJacobian>(1); - testAlphaMinusJacobian>(2); - testAlphaMinusJacobian>(3); - testAlphaMinusJacobian>(4); + testAlphaMinusJacobian>(1); + testAlphaMinusJacobian>(2); + testAlphaMinusJacobian>(3); + testAlphaMinusJacobian>(4); } template @@ -168,29 +167,19 @@ using Group8VectorMatrix = micm::VectorMatrix; template using Group10VectorMatrix = micm::VectorMatrix; -template -using Group1SparseVectorMatrix = micm::SparseMatrix>; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; -template -using Group3SparseVectorMatrix = micm::SparseMatrix>; -template -using Group4SparseVectorMatrix = micm::SparseMatrix>; -template -using Group8SparseVectorMatrix = micm::SparseMatrix>; -template -using Group10SparseVectorMatrix = micm::SparseMatrix>; +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; +using Group8SparseVectorMatrix = micm::SparseMatrix>; +using Group10SparseVectorMatrix = micm::SparseMatrix>; TEST(RosenbrockSolver, DenseAlphaMinusJacobian) { - testAlphaMinusJacobian>( - 1); - testAlphaMinusJacobian>( - 4); - testAlphaMinusJacobian>( - 3); - testAlphaMinusJacobian>( - 2); + testAlphaMinusJacobian>(1); + testAlphaMinusJacobian>(4); + testAlphaMinusJacobian>(3); + testAlphaMinusJacobian>(2); } TEST(RosenbrockSolver, CanSetTolerances) @@ -252,42 +241,15 @@ TEST(RosenbrockSolver, CanOverrideTolerancesWithParameters) TEST(RosenbrockSolver, NormalizedError) { // Exact fits - testNormalizedErrorDiff< - Group1VectorMatrix, - Group1SparseVectorMatrix, - micm::LinearSolver>(1); - testNormalizedErrorDiff< - Group2VectorMatrix, - Group2SparseVectorMatrix, - micm::LinearSolver>(2); - testNormalizedErrorDiff< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::LinearSolver>(3); - testNormalizedErrorDiff< - Group4VectorMatrix, - Group4SparseVectorMatrix, - micm::LinearSolver>(4); + testNormalizedErrorDiff>(1); + testNormalizedErrorDiff>(2); + testNormalizedErrorDiff>(3); + testNormalizedErrorDiff>(4); // Inexact fits - testNormalizedErrorDiff< - Group2VectorMatrix, - Group2SparseVectorMatrix, - micm::LinearSolver>(1); - testNormalizedErrorDiff< - Group3VectorMatrix, - Group3SparseVectorMatrix, - micm::LinearSolver>(2); - testNormalizedErrorDiff< - Group4VectorMatrix, - Group4SparseVectorMatrix, - micm::LinearSolver>(3); - testNormalizedErrorDiff< - Group8VectorMatrix, - Group8SparseVectorMatrix, - micm::LinearSolver>(5); - testNormalizedErrorDiff< - Group10VectorMatrix, - Group10SparseVectorMatrix, - micm::LinearSolver>(3); + testNormalizedErrorDiff>(1); + testNormalizedErrorDiff>(2); + testNormalizedErrorDiff>(3); + testNormalizedErrorDiff>(5); + testNormalizedErrorDiff>(3); } \ No newline at end of file diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp new file mode 100644 index 000000000..28f60b141 --- /dev/null +++ b/test/unit/solver/test_solver_builder.cpp @@ -0,0 +1,96 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace +{ + auto a = micm::Species("A"); + auto b = micm::Species("B"); + auto c = micm::Species("C"); + + micm::Phase gas_phase{ std::vector{ a, b, c } }; + + micm::Process r1 = micm::Process::Create() + .SetReactants({ a }) + .SetProducts({ micm::Yields(b, 1) }) + .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 2.15e-11, .B_ = 0, .C_ = 110 })) + .SetPhase(gas_phase); + + micm::Process r2 = micm::Process::Create() + .SetReactants({ b }) + .SetProducts({ micm::Yields(c, 1) }) + .SetRateConstant(micm::ArrheniusRateConstant( + micm::ArrheniusRateConstantParameters{ .A_ = 3.3e-11, .B_ = 0, .C_ = 55 })) + .SetPhase(gas_phase); + micm::System the_system = micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }); + std::vector reactions = { r1, r2 }; +} // namespace + +TEST(SolverBuilder, CanBuildBackwardEuler) +{ + auto backward_euler = micm::CpuSolverBuilder, micm::SparseMatrix>() + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(1) + .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) + .SetSolverParameters(micm::BackwardEulerSolverParameters{}) + .Build(); +} + +TEST(SolverBuilder, CanBuildRosenbrock) +{ + // auto rosenbrock = micm::CpuSolverBuilder, SparseMatrix>() + // .SetSystem(the_system) + // .SetReactions(reactions) + // .SetNumberOfGridCells(1) + // .SolverParameters(micm::ThreeStageRosenbockSolverParameters{}) + // .Build(); + + // auto rosenbrock_vector = micm::CpuSolverBuilder, SparseVectorMatrix>() + // .SetSystem(the_system) + // .SetReactions(reactions) + // .SetNumberOfGridCells(1) + // .SolverParameters(micm::ThreeStageRosenbockSolverParameters{}) + // .Build(); +} + +TEST(SolverBuilder, CanBuildJitRosenbrock) +{ + // auto jit_rosenbrock = micm::JitSolverBuilder() + // .SetSystem(the_system) + // .SetReactions(reactions) + // .SetNumberOfGridCells(1) + // .SolverParameters(micm::ThreeStageRosenbockSolverParameters{}) + // .Build(); +} + +TEST(SolverBuilder, CanBuildCudaSolvers) +{ + // auto cuda_backward_euler = micm::CudaSolverBuilder() + // .SetSystem(the_system) + // .SetReactions(reactions) + // .SetNumberOfGridCells(1) + // .SolverParameters(micm::BackwardEulerSolverParameters{}) + // .Build(); + + // auto cuda_rosenbrock = micm::CudaSolverBuilder() + // .SetSystem(the_system) + // .SetReactions(reactions) + // .SetNumberOfGridCells(1) + // .SolverParameters(micm::ThreeStageRosenbockSolverParameters{}) + // .Build(); +} \ No newline at end of file diff --git a/test/valgrind.supp b/test/valgrind.supp index e45889e52..99dc72fe7 100644 --- a/test/valgrind.supp +++ b/test/valgrind.supp @@ -10,4 +10,14 @@ fun:_ZN4llvm10AsmPrinter28emitPatchableFunctionEntriesEv fun:_ZN4llvm10AsmPrinter16emitFunctionBodyEv ... +} + +{ + + Memcheck:Leak + match-leak-kinds: possible + fun:calloc + ... + fun:GOMP_parallel + ... } \ No newline at end of file