diff --git a/CMakeLists.txt b/CMakeLists.txt index 8f19d7a3..175f4b8b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,6 +30,8 @@ option(MUSICA_BUILD_DOCS "Build the documentation" OFF) option(MUSICA_ENABLE_MICM "Enable MICM" ON) option(MUSICA_ENABLE_TUVX "Enable TUV-x" ON) +set(MUSICA_SET_MICM_VECTOR_MATRIX_SIZE "1" CACHE STRING "Set MICM vector-ordered matrix dimension") + cmake_dependent_option( MUSICA_ENABLE_PYTHON_LIBRARY "Adds pybind11, a lightweight header-only library that exposes C++ types in Python and vice versa" OFF "MUSICA_BUILD_C_CXX_INTERFACE" OFF) diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake index 9b44a537..df84a956 100644 --- a/cmake/dependencies.cmake +++ b/cmake/dependencies.cmake @@ -60,6 +60,7 @@ if (MUSICA_ENABLE_MICM AND MUSICA_BUILD_C_CXX_INTERFACE) ) set(MICM_ENABLE_TESTS OFF) set(MICM_ENABLE_EXAMPLES OFF) + set(MICM_DEFAULT_VECTOR_MATRIX_SIZE ${MUSICA_SET_MICM_VECTOR_MATRIX_SIZE}) FetchContent_MakeAvailable(micm) endif() diff --git a/cmake/test_util.cmake b/cmake/test_util.cmake index 074f01d5..39cafea1 100644 --- a/cmake/test_util.cmake +++ b/cmake/test_util.cmake @@ -43,7 +43,7 @@ function(create_standard_test_cxx) add_executable(test_${TEST_NAME} ${TEST_SOURCES}) target_link_libraries(test_${TEST_NAME} PUBLIC musica::musica GTest::gtest_main) if(MUSICA_ENABLE_OPENMP) - target_link_libraries(test_${TEST_NAME} PUBLIC OpenMP::OpenMP_CXX OpenMP::OpenMP_Fortran) + target_link_libraries(test_${TEST_NAME} PUBLIC OpenMP::OpenMP_CXX) endif() if(NOT DEFINED TEST_WORKING_DIRECTORY) set(TEST_WORKING_DIRECTORY "${CMAKE_BINARY_DIR}") diff --git a/docker/Dockerfile.fortran-gcc b/docker/Dockerfile.fortran-gcc index 7e6410c9..d21acb7e 100644 --- a/docker/Dockerfile.fortran-gcc +++ b/docker/Dockerfile.fortran-gcc @@ -1,6 +1,6 @@ FROM fedora:35 -ARG BUILD_TYPE=release +ARG BUILD_TYPE=Release RUN dnf -y update \ && dnf -y install \ diff --git a/docker/Dockerfile.fortran-intel b/docker/Dockerfile.fortran-intel index 3e0dac02..34bb0eb7 100644 --- a/docker/Dockerfile.fortran-intel +++ b/docker/Dockerfile.fortran-intel @@ -51,7 +51,7 @@ COPY . musica RUN cd musica \ && cmake -S . \ -B build \ - -D CMAKE_BUILD_TYPE=Release \ + -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \ && cd build \ && make install -j @@ -69,6 +69,4 @@ RUN cd musica/fortran/test/fetch_content_integration \ -D MUSICA_ENABLE_MEMCHECK=ON \ && make -j -WORKDIR musica/fortran/test/fetch_content_integration/build -RUN cp -r /musica/build/_deps/tuvx-src/examples/ . -RUN cp -r /musica/build/_deps/tuvx-src/data/ . \ No newline at end of file +WORKDIR musica/fortran/test/fetch_content_integration/build \ No newline at end of file diff --git a/docker/Dockerfile.memcheck b/docker/Dockerfile.memcheck index 25c61637..cd943615 100644 --- a/docker/Dockerfile.memcheck +++ b/docker/Dockerfile.memcheck @@ -1,5 +1,7 @@ FROM fedora:35 +ARG BUILD_TYPE=Debug + RUN dnf -y update \ && dnf -y install \ cmake \ @@ -22,7 +24,7 @@ COPY . musica RUN cd musica \ && cmake -S . \ -B build \ - -D CMAKE_BUILD_TYPE=Debug \ + -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \ -D MUSICA_ENABLE_MEMCHECK=ON \ && cd build \ && make install -j 8 diff --git a/docker/Dockerfile.mpi b/docker/Dockerfile.mpi index 91acf76e..fa85b48f 100644 --- a/docker/Dockerfile.mpi +++ b/docker/Dockerfile.mpi @@ -1,5 +1,7 @@ FROM fedora:35 +ARG BUILD_TYPE=Debug + RUN dnf -y update \ && dnf install -y sudo \ && adduser test_user \ @@ -38,9 +40,9 @@ RUN sudo chown -R test_user.test_user musica RUN cd musica \ && cmake -S . \ -B build \ - -D CMAKE_BUILD_TYPE=debug \ - -D ENABLE_TESTS=ON \ - -D ENABLE_MPI=ON \ + -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \ + -D MUSICA_ENABLE_TESTS=ON \ + -D MUSICA_ENABLE_MPI=ON \ -D CMAKE_Fortran_COMPILER=/usr/lib64/openmpi/bin/mpif90 \ -D CMAKE_C_COMPILER=/usr/lib64/openmpi/bin/mpicc \ -D CMAKE_CXX_COMPILER=/usr/lib64/openmpi/bin/mpicxx \ diff --git a/docker/Dockerfile.mpi_openmp b/docker/Dockerfile.mpi_openmp index a9147545..f48a4605 100644 --- a/docker/Dockerfile.mpi_openmp +++ b/docker/Dockerfile.mpi_openmp @@ -1,5 +1,7 @@ FROM fedora:35 +ARG BUILD_TYPE=Debug + RUN dnf -y update \ && dnf install -y sudo \ && adduser test_user \ @@ -38,10 +40,10 @@ RUN sudo chown -R test_user.test_user musica RUN cd musica \ && cmake -S . \ -B build \ - -D CMAKE_BUILD_TYPE=debug \ - -D ENABLE_MPI=ON \ - -D ENABLE_OPENMP=ON \ - -D ENABLE_TESTS=ON \ + -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \ + -D MUSICA_ENABLE_MPI=ON \ + -D MUSICA_ENABLE_OPENMP=ON \ + -D MUSICA_ENABLE_TESTS=ON \ -D CMAKE_Fortran_COMPILER=/usr/lib64/openmpi/bin/mpif90 \ -D CMAKE_C_COMPILER=/usr/lib64/openmpi/bin/mpicc \ -D CMAKE_CXX_COMPILER=/usr/lib64/openmpi/bin/mpicxx \ diff --git a/docker/Dockerfile.openmp b/docker/Dockerfile.openmp index 7ad2c7af..efc433b8 100644 --- a/docker/Dockerfile.openmp +++ b/docker/Dockerfile.openmp @@ -1,5 +1,7 @@ FROM fedora:35 +ARG BUILD_TYPE=Debug + RUN dnf -y update \ && dnf install -y sudo \ && adduser test_user \ @@ -37,9 +39,9 @@ RUN sudo chown -R test_user.test_user musica RUN cd musica \ && cmake -S . \ -B build \ - -D CMAKE_BUILD_TYPE=debug \ - -D ENABLE_OPENMP:BOOL=TRUE \ - -D ENABLE_TESTS=ON \ + -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \ + -D MUSICA_ENABLE_OPENMP=ON \ + -D MUSICA_ENABLE_TESTS=ON \ -D CMAKE_Fortran_COMPILER=/usr/lib64/openmpi/bin/mpif90 \ -D CMAKE_C_COMPILER=/usr/lib64/openmpi/bin/mpicc \ -D CMAKE_CXX_COMPILER=/usr/lib64/openmpi/bin/mpic++ \ diff --git a/docker/Dockerfile.python b/docker/Dockerfile.python index e7c27299..02788a0f 100644 --- a/docker/Dockerfile.python +++ b/docker/Dockerfile.python @@ -1,5 +1,7 @@ FROM fedora:latest +ARG BUILD_TYPE=Release + RUN dnf -y update \ && dnf -y install \ cmake \ @@ -45,7 +47,7 @@ RUN pip install -r requirements.txt RUN cd musica \ && cmake -S . \ -B build \ - -D CMAKE_BUILD_TYPE=Release \ + -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \ -D MUSICA_ENABLE_PYTHON_LIBRARY=ON \ -D MUSICA_ENABLE_TUVX=OFF \ && cd build \ diff --git a/fortran/micm.F90 b/fortran/micm.F90 index 854a4db8..db1ab198 100644 --- a/fortran/micm.F90 +++ b/fortran/micm.F90 @@ -10,6 +10,7 @@ module musica_micm implicit none public :: micm_t, solver_stats_t, get_micm_version + public :: Rosenbrock, RosenbrockStandardOrder private !> Wrapper for c solver stats @@ -25,13 +26,22 @@ module musica_micm real(c_double) :: final_time_ = 0._c_double end type solver_stats_t_c + ! We could use Fortran 2023 enum type feature if Fortran 2023 is supported + ! https://fortran-lang.discourse.group/t/enumerator-type-in-bind-c-derived-type-best-practice/5947/2 + enum, bind(c) + enumerator :: Rosenbrock = 1 + enumerator :: RosenbrockStandardOrder = 2 + end enum + interface - function create_micm_c(config_path, error) bind(C, name="CreateMicm") + function create_micm_c(config_path, solver_type, num_grid_cells, error) bind(C, name="CreateMicm") use musica_util, only: error_t_c import c_ptr, c_int, c_char - character(kind=c_char), intent(in) :: config_path(*) - type(error_t_c), intent(inout) :: error - type(c_ptr) :: create_micm_c + character(kind=c_char), intent(in) :: config_path(*) + integer(kind=c_int), value, intent(in) :: solver_type + integer(kind=c_int), value, intent(in) :: num_grid_cells + type(error_t_c), intent(inout) :: error + type(c_ptr) :: create_micm_c end function create_micm_c subroutine delete_micm_c(micm, error) bind(C, name="DeleteMicm") @@ -182,10 +192,12 @@ function get_micm_version() result(value) value = string_t(string_c) end function get_micm_version - function constructor(config_path, error) result( this ) + function constructor(config_path, solver_type, num_grid_cells, error) result( this ) use musica_util, only: error_t_c, error_t, copy_mappings type(micm_t), pointer :: this character(len=*), intent(in) :: config_path + integer(c_int), intent(in) :: solver_type + integer(c_int), intent(in) :: num_grid_cells type(error_t), intent(inout) :: error character(len=1, kind=c_char) :: c_config_path(len_trim(config_path)+1) integer :: n, i @@ -201,7 +213,7 @@ function constructor(config_path, error) result( this ) end do c_config_path(n+1) = c_null_char - this%ptr = create_micm_c(c_config_path, error_c) + this%ptr = create_micm_c(c_config_path, solver_type, num_grid_cells, error_c) error = error_t(error_c) if (.not. error%is_success()) then deallocate(this) diff --git a/fortran/test/fetch_content_integration/test_micm_api.F90 b/fortran/test/fetch_content_integration/test_micm_api.F90 index c98eca9d..24d9d31c 100644 --- a/fortran/test/fetch_content_integration/test_micm_api.F90 +++ b/fortran/test/fetch_content_integration/test_micm_api.F90 @@ -6,6 +6,7 @@ program test_micm_api use, intrinsic :: iso_c_binding use, intrinsic :: ieee_arithmetic use musica_micm, only: micm_t, solver_stats_t, get_micm_version + use musica_micm, only: Rosenbrock, RosenbrockStandardOrder use musica_util, only: assert, error_t, mapping_t, string_t #include "micm/util/error.hpp" @@ -32,6 +33,8 @@ subroutine test_api() real(c_double), dimension(5) :: concentrations real(c_double), dimension(3) :: user_defined_reaction_rates character(len=256) :: config_path + integer(c_int) :: solver_type + integer(c_int) :: num_grid_cells character(len=:), allocatable :: string_value real(c_double) :: double_value integer(c_int) :: int_value @@ -42,13 +45,15 @@ subroutine test_api() real(c_double), parameter :: GAS_CONSTANT = 8.31446261815324_c_double ! J mol-1 K-1 integer :: i + config_path = "configs/chapman" + solver_type = Rosenbrock + num_grid_cells = 1 time_step = 200 temperature = 272.5 pressure = 101253.4 air_density = pressure / ( GAS_CONSTANT * temperature ) num_concentrations = 5 concentrations = (/ 0.75, 0.4, 0.8, 0.01, 0.02 /) - config_path = "configs/chapman" num_user_defined_reaction_rates = 3 user_defined_reaction_rates = (/ 0.1, 0.2, 0.3 /) @@ -56,7 +61,7 @@ subroutine test_api() print *, "[test micm fort api] MICM version ", micm_version%get_char_array() write(*,*) "[test micm fort api] Creating MICM solver..." - micm => micm_t(config_path, error) + micm => micm_t(config_path, solver_type, num_grid_cells, error) ASSERT( error%is_success() ) do i = 1, size( micm%species_ordering ) @@ -116,7 +121,7 @@ subroutine test_api() ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, \ MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) ) deallocate( micm ) - micm => micm_t( "configs/invalid", error ) + micm => micm_t( "configs/invalid", solver_type, num_grid_cells, error ) ASSERT( error%is_error( MICM_ERROR_CATEGORY_CONFIGURATION, \ MICM_CONFIGURATION_ERROR_CODE_INVALID_FILE_PATH ) ) ASSERT( .not. associated( micm ) ) diff --git a/fortran/test/fetch_content_integration/test_micm_box_model.F90 b/fortran/test/fetch_content_integration/test_micm_box_model.F90 index dfdc5aab..b80a4af1 100644 --- a/fortran/test/fetch_content_integration/test_micm_box_model.F90 +++ b/fortran/test/fetch_content_integration/test_micm_box_model.F90 @@ -5,6 +5,7 @@ program test_micm_box_model use musica_util, only: error_t, string_t, mapping_t use musica_micm, only: micm_t, solver_stats_t + use musica_micm, only: Rosenbrock, RosenbrockStandardOrder implicit none @@ -15,6 +16,8 @@ program test_micm_box_model subroutine box_model() character(len=256) :: config_path + integer(c_int) :: solver_type + integer(c_int) :: num_grid_cells real(c_double), parameter :: GAS_CONSTANT = 8.31446261815324_c_double ! J mol-1 K-1 @@ -38,6 +41,8 @@ subroutine box_model() integer :: i config_path = "configs/analytical" + solver_type = RosenbrockStandardOrder + num_grid_cells = 1 time_step = 200 temperature = 273.0 @@ -47,7 +52,7 @@ subroutine box_model() concentrations = (/ 1.0, 1.0, 1.0 /) write(*,*) "Creating MICM solver..." - micm => micm_t(config_path, error) + micm => micm_t(config_path, solver_type, num_grid_cells, error) do i = 1, size( micm%species_ordering ) associate(the_mapping => micm%species_ordering(i)) diff --git a/include/musica/micm.hpp b/include/musica/micm.hpp index c4b10075..a7356f6e 100644 --- a/include/musica/micm.hpp +++ b/include/musica/micm.hpp @@ -12,13 +12,21 @@ #include #include #include +#include #include +#include +#include #include +#include #include #include #include +#ifndef MICM_VECTOR_MATRIX_SIZE + #define MICM_VECTOR_MATRIX_SIZE 1 +#endif + namespace musica { @@ -28,6 +36,12 @@ namespace musica extern "C" { #endif + /// @brief Types of MICM solver + enum MICMSolver + { + Rosenbrock = 1, // Vector-ordered Rosenbrock solver + RosenbrockStandardOrder, // Standard-ordered Rosenbrock solver + }; struct SolverResultStats { @@ -87,7 +101,12 @@ namespace musica } }; - MICM *CreateMicm(const char *config_path, Error *error); + /// @brief Create a MICM object by specifying solver type to use + /// @param config_path Path to configuration file or directory containing configuration file + /// @param solver_type Type of MICMSolver + /// @param num_grid_cells Number of grid cells + /// @param error Error struct to indicate success or failure + MICM *CreateMicm(const char *config_path, MICMSolver solver_type, int num_grid_cells, Error *error); void DeleteMicm(const MICM *micm, Error *error); void MicmSolve( MICM *micm, @@ -116,12 +135,18 @@ namespace musica class MICM { public: - /// @brief Create a solver by reading and parsing configuration file + /// @brief Create a Rosenbrock solver of vector-ordered matrix type by reading and parsing configuration file + /// @param config_path Path to configuration file or directory containing configuration file + /// @param error Error struct to indicate success or failure + void CreateRosenbrock(const std::string &config_path, Error *error); + + /// @brief Create a Rosenbrock solver of standard-ordered matrix type by reading and parsing configuration file /// @param config_path Path to configuration file or directory containing configuration file /// @param error Error struct to indicate success or failure - void Create(const std::string &config_path, Error *error); + void CreateRosenbrockStandardOrder(const std::string &config_path, Error *error); /// @brief Solve the system + /// @param solver Pointer to solver /// @param time_step Time [s] to advance the state by /// @param temperature Temperature [K] /// @param pressure Pressure [Pa] @@ -132,6 +157,7 @@ namespace musica /// @param custom_rate_parameters Array of custom rate parameters /// @param error Error struct to indicate success or failure void Solve( + auto &solver, double time_step, double temperature, double pressure, @@ -144,15 +170,34 @@ namespace musica SolverResultStats *solver_stats, Error *error); + /// @brief Set solver type + /// @param MICMSolver Type of MICMSolver + void SetSolverType(MICMSolver solver_type) + { + solver_type_ = solver_type; + } + + /// @brief Set number of grid cells + /// @param num_grid_cells Number of grid cells + void SetNumGridCells(int num_grid_cells) + { + num_grid_cells_ = num_grid_cells; + } + /// @brief Get the ordering of species + /// @param solver Pointer to solver /// @param error Error struct to indicate success or failure /// @return Map of species names to their indices - std::map GetSpeciesOrdering(Error *error); + // std::map GetSpeciesOrdering(auto &solver, Error *error); + template + std::map GetSpeciesOrdering(T &solver, Error *error); /// @brief Get the ordering of user-defined reaction rates + /// @param solver Pointer to solver /// @param error Error struct to indicate success or failure /// @return Map of reaction rate names to their indices - std::map GetUserDefinedReactionRatesOrdering(Error *error); + template + std::map GetUserDefinedReactionRatesOrdering(T &solver, Error *error); /// @brief Get a property for a chemical species /// @param species_name Name of the species @@ -162,20 +207,66 @@ namespace musica template T GetSpeciesProperty(const std::string &species_name, const std::string &property_name, Error *error); - static constexpr std::size_t NUM_GRID_CELLS = 1; + public: + MICMSolver solver_type_; - private: - using DenseMatrixPolicy = micm::Matrix; - using SparseMatrixPolicy = micm::SparseMatrix; - using SolverPolicy = typename micm::RosenbrockSolverParameters:: - template SolverType>; - using Rosenbrock = micm::Solver>; + /// @brief Vector-ordered Rosenbrock solver type + using DenseMatrixVector = micm::VectorMatrix; + using SparseMatrixVector = micm::SparseMatrix>; + using RosenbrockVectorType = typename micm::RosenbrockSolverParameters:: + template SolverType>; + using Rosenbrock = micm::Solver>; + std::unique_ptr rosenbrock_; - std::unique_ptr solver_; + /// @brief Standard-ordered Rosenbrock solver type + using DenseMatrixStandard = micm::Matrix; + using SparseMatrixStandard = micm::SparseMatrix; + using RosenbrockStandardType = typename micm::RosenbrockSolverParameters:: + template SolverType>; + using RosenbrockStandard = micm::Solver>; + std::unique_ptr rosenbrock_standard_; + private: + int num_grid_cells_; std::unique_ptr solver_parameters_; }; + template + inline std::map MICM::GetSpeciesOrdering(T &solver, Error *error) + { + try + { + micm::State state = solver->GetState(); + DeleteError(error); + *error = NoError(); + return state.variable_map_; + } + catch (const std::system_error &e) + { + DeleteError(error); + *error = ToError(e); + return std::map(); + } + } + + template + inline std::map MICM::GetUserDefinedReactionRatesOrdering(T &solver, Error *error) + { + try + { + micm::State state = solver->GetState(); + DeleteError(error); + *error = NoError(); + return state.custom_rate_parameter_map_; + } + catch (const std::system_error &e) + { + DeleteError(error); + *error = ToError(e); + return std::map(); + } + } + template inline T MICM::GetSpeciesProperty(const std::string &species_name, const std::string &property_name, Error *error) { @@ -201,4 +292,4 @@ namespace musica *error = ToError(MUSICA_ERROR_CATEGORY, MUSICA_ERROR_CODE_SPECIES_NOT_FOUND, msg.c_str()); return T(); } -} // namespace musica +} // namespace musica \ No newline at end of file diff --git a/include/musica/util.hpp b/include/musica/util.hpp index e6f3f0c6..aaca47b0 100644 --- a/include/musica/util.hpp +++ b/include/musica/util.hpp @@ -4,8 +4,9 @@ #include -#define MUSICA_ERROR_CATEGORY "MUSICA Error" -#define MUSICA_ERROR_CODE_SPECIES_NOT_FOUND 1 +#define MUSICA_ERROR_CATEGORY "MUSICA Error" +#define MUSICA_ERROR_CODE_SPECIES_NOT_FOUND 1 +#define MUSICA_ERROR_CODE_SOLVER_TYPE_NOT_FOUND 2 #ifdef __cplusplus #include diff --git a/python/test/test_analytical.py b/python/test/test_analytical.py index 6f973847..d446447f 100644 --- a/python/test/test_analytical.py +++ b/python/test/test_analytical.py @@ -5,14 +5,14 @@ class TestAnalyticalSimulation(unittest.TestCase): def test_simulation(self): + num_grid_cells = 1 time_step = 200.0 temperature = 272.5 pressure = 101253.3 GAS_CONSTANT = 8.31446261815324 air_density = pressure / (GAS_CONSTANT * temperature) - solver = musica.create_solver("configs/analytical") - + solver = musica.create_solver("configs/analytical", musica.micmsolver.rosenbrock, num_grid_cells) rates = musica.user_defined_reaction_rates(solver) ordering = musica.species_ordering(solver) diff --git a/python/test/test_chapman.py b/python/test/test_chapman.py index 90e93acf..fdfd2745 100644 --- a/python/test/test_chapman.py +++ b/python/test/test_chapman.py @@ -4,6 +4,7 @@ class TestChapman(unittest.TestCase): def test_micm_solve(self): + num_grid_cells = 1 time_step = 200.0 temperature = 272.5 pressure = 101253.3 @@ -11,7 +12,7 @@ def test_micm_solve(self): air_density = pressure / (GAS_CONSTANT * temperature) concentrations = [0.75, 0.4, 0.8, 0.01, 0.02] - solver = musica.create_solver("configs/chapman") + solver = musica.create_solver("configs/chapman", musica.micmsolver.rosenbrock, num_grid_cells) rate_constant_ordering = musica.user_defined_reaction_rates(solver) ordering = musica.species_ordering(solver) diff --git a/python/wrapper.cpp b/python/wrapper.cpp index 7310baf9..8f828f19 100644 --- a/python/wrapper.cpp +++ b/python/wrapper.cpp @@ -10,18 +10,20 @@ namespace py = pybind11; // Wraps micm.cpp PYBIND11_MODULE(musica, m) { - py::class_(m, "MICM") + py::class_(m, "micm") .def(py::init<>()) - .def("create", &musica::MICM::Create) - .def("solve", &musica::MICM::Solve) .def("__del__", [](musica::MICM &micm) {}); + py::enum_(m, "micmsolver") + .value("rosenbrock", musica::MICMSolver::Rosenbrock) + .value("rosenbrock_standard_order", musica::MICMSolver::RosenbrockStandardOrder); + m.def( "create_solver", - [](const char *config_path) + [](const char *config_path, musica::MICMSolver solver_type, int num_grid_cells) { musica::Error error; - musica::MICM *micm = musica::CreateMicm(config_path, &error); + musica::MICM *micm = musica::CreateMicm(config_path, solver_type, num_grid_cells, &error); if (!musica::IsSuccess(error)) { std::string message = "Error creating solver: " + std::string(error.message_.value_); @@ -90,7 +92,18 @@ PYBIND11_MODULE(musica, m) [](musica::MICM *micm) { musica::Error error; - return micm->GetSpeciesOrdering(&error); + std::map map; + + if (micm->solver_type_ == musica::MICMSolver::Rosenbrock) + { + map = micm->GetSpeciesOrdering(micm->rosenbrock_, &error); + } + else if (micm->solver_type_ == musica::MICMSolver::RosenbrockStandardOrder) + { + map = micm->GetSpeciesOrdering(micm->rosenbrock_standard_, &error); + } + + return map; }, "Return map of get_species_ordering rates"); @@ -99,7 +112,18 @@ PYBIND11_MODULE(musica, m) [](musica::MICM *micm) { musica::Error error; - return micm->GetUserDefinedReactionRatesOrdering(&error); + std::map map; + + if (micm->solver_type_ == musica::MICMSolver::Rosenbrock) + { + map = micm->GetUserDefinedReactionRatesOrdering(micm->rosenbrock_, &error); + } + else if (micm->solver_type_ == musica::MICMSolver::RosenbrockStandardOrder) + { + map = micm->GetUserDefinedReactionRatesOrdering(micm->rosenbrock_standard_, &error); + } + + return map; }, "Return map of reaction rates"); } \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e8c7a849..54a5f60f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -17,6 +17,8 @@ message (STATUS "CMake build configuration for ${PROJECT_NAME} (${CMAKE_BUILD_TY add_library(musica) add_library(musica::musica ALIAS musica) +target_compile_definitions(musica PUBLIC MICM_VECTOR_MATRIX_SIZE=${MUSICA_SET_MICM_VECTOR_MATRIX_SIZE}) + # set the c++ standard for musica target_compile_features(musica PUBLIC cxx_std_20) diff --git a/src/micm/micm.cpp b/src/micm/micm.cpp index 7d03038e..dd5b804a 100644 --- a/src/micm/micm.cpp +++ b/src/micm/micm.cpp @@ -5,6 +5,7 @@ // multi-component reactive transport model. It also includes functions for // creating and deleting MICM instances, creating solvers, and solving the model. #include +#include #include #include @@ -14,20 +15,41 @@ #include #include #include -#include +#include +#include namespace musica { - MICM *CreateMicm(const char *config_path, Error *error) + + MICM *CreateMicm(const char *config_path, MICMSolver solver_type, int num_grid_cells, Error *error) { DeleteError(error); MICM *micm = new MICM(); - micm->Create(std::string(config_path), error); + micm->SetNumGridCells(num_grid_cells); + + if (solver_type == MICMSolver::Rosenbrock) + { + micm->SetSolverType(MICMSolver::Rosenbrock); + micm->CreateRosenbrock(std::string(config_path), error); + } + else if (solver_type == MICMSolver::RosenbrockStandardOrder) + { + micm->SetSolverType(MICMSolver::RosenbrockStandardOrder); + micm->CreateRosenbrockStandardOrder(std::string(config_path), error); + } + else + { + std::string msg = "Solver type '" + std::to_string(solver_type) + "' not found"; + *error = ToError(MUSICA_ERROR_CATEGORY, MUSICA_ERROR_CODE_SOLVER_TYPE_NOT_FOUND, msg.c_str()); + delete micm; + return nullptr; + } if (!IsSuccess(*error)) { delete micm; return nullptr; } + return micm; } @@ -65,19 +87,40 @@ namespace musica Error *error) { DeleteError(error); - micm->Solve( - time_step, - temperature, - pressure, - air_density, - num_concentrations, - concentrations, - num_custom_rate_parameters, - custom_rate_parameters, - solver_state, - solver_stats, - error); - } + + if (micm->solver_type_ == MICMSolver::Rosenbrock) + { + micm->Solve( + micm->rosenbrock_, + time_step, + temperature, + pressure, + air_density, + num_concentrations, + concentrations, + num_custom_rate_parameters, + custom_rate_parameters, + solver_state, + solver_stats, + error); + } + else if (micm->solver_type_ == MICMSolver::RosenbrockStandardOrder) + { + micm->Solve( + micm->rosenbrock_standard_, + time_step, + temperature, + pressure, + air_density, + num_concentrations, + concentrations, + num_custom_rate_parameters, + custom_rate_parameters, + solver_state, + solver_stats, + error); + } + }; String MicmVersion() { @@ -87,11 +130,22 @@ namespace musica Mapping *GetSpeciesOrdering(MICM *micm, std::size_t *array_size, Error *error) { DeleteError(error); - auto map = micm->GetSpeciesOrdering(error); + + std::map map; + + if (micm->solver_type_ == MICMSolver::Rosenbrock) + { + map = micm->GetSpeciesOrdering(micm->rosenbrock_, error); + } + else if (micm->solver_type_ == MICMSolver::RosenbrockStandardOrder) + { + map = micm->GetSpeciesOrdering(micm->rosenbrock_standard_, error); + } if (!IsSuccess(*error)) { return nullptr; } + Mapping *species_ordering = new Mapping[map.size()]; // Copy data from the map to the array of structs @@ -110,11 +164,22 @@ namespace musica Mapping *GetUserDefinedReactionRatesOrdering(MICM *micm, std::size_t *array_size, Error *error) { DeleteError(error); - auto map = micm->GetUserDefinedReactionRatesOrdering(error); + + std::map map; + + if (micm->solver_type_ == MICMSolver::Rosenbrock) + { + map = micm->GetUserDefinedReactionRatesOrdering(micm->rosenbrock_, error); + } + else if (micm->solver_type_ == MICMSolver::RosenbrockStandardOrder) + { + map = micm->GetUserDefinedReactionRatesOrdering(micm->rosenbrock_standard_, error); + } if (!IsSuccess(*error)) { return nullptr; } + Mapping *reactionRates = new Mapping[map.size()]; // Copy data from the map to the array of structs @@ -167,7 +232,7 @@ namespace musica return micm->GetSpeciesProperty(species_name_str, property_name_str, error); } - void MICM::Create(const std::string &config_path, Error *error) + void MICM::CreateRosenbrock(const std::string &config_path, Error *error) { try { @@ -175,13 +240,47 @@ namespace musica solver_config.ReadAndParse(std::filesystem::path(config_path)); solver_parameters_ = std::make_unique(solver_config.GetSolverParams()); - solver_ = std::make_unique(micm::CpuSolverBuilder( - micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) - .SetSystem(solver_parameters_->system_) - .SetReactions(solver_parameters_->processes_) - .SetNumberOfGridCells(NUM_GRID_CELLS) - .SetIgnoreUnusedSpecies(true) - .Build()); + rosenbrock_ = std::make_unique( + micm::SolverBuilder< + micm::RosenbrockSolverParameters, + micm::VectorMatrix, + micm::SparseMatrix>, + micm::ProcessSet, + micm::LinearSolver< + micm::SparseMatrix>, + micm::LuDecomposition>>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(solver_parameters_->system_) + .SetReactions(solver_parameters_->processes_) + .SetNumberOfGridCells(num_grid_cells_) + .SetIgnoreUnusedSpecies(true) + .Build()); + + DeleteError(error); + *error = NoError(); + } + catch (const std::system_error &e) + { + DeleteError(error); + *error = ToError(e); + } + } + + void MICM::CreateRosenbrockStandardOrder(const std::string &config_path, Error *error) + { + try + { + micm::SolverConfig<> solver_config; + solver_config.ReadAndParse(std::filesystem::path(config_path)); + solver_parameters_ = std::make_unique(solver_config.GetSolverParams()); + + rosenbrock_standard_ = + std::make_unique(micm::CpuSolverBuilder( + micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(solver_parameters_->system_) + .SetReactions(solver_parameters_->processes_) + .SetNumberOfGridCells(num_grid_cells_) + .SetIgnoreUnusedSpecies(true) + .Build()); DeleteError(error); *error = NoError(); @@ -194,6 +293,7 @@ namespace musica } void MICM::Solve( + auto &solver, double time_step, double temperature, double pressure, @@ -208,13 +308,13 @@ namespace musica { try { - micm::State state = solver_->GetState(); + micm::State state = solver->GetState(); - for (std::size_t i{}; i < NUM_GRID_CELLS; i++) + for (int cell{}; cell < num_grid_cells_; cell++) { - state.conditions_[i].temperature_ = temperature; - state.conditions_[i].pressure_ = pressure; - state.conditions_[i].air_density_ = air_density; + state.conditions_[cell].temperature_ = temperature; + state.conditions_[cell].pressure_ = pressure; + state.conditions_[cell].air_density_ = air_density; } state.variables_.AsVector().assign(concentrations, concentrations + num_concentrations); @@ -222,8 +322,8 @@ namespace musica state.custom_rate_parameters_.AsVector().assign( custom_rate_parameters, custom_rate_parameters + num_custom_rate_parameters); - solver_->CalculateRateConstants(state); - auto result = solver_->Solve(time_step, state); + solver->CalculateRateConstants(state); + auto result = solver->Solve(time_step, state); *solver_state = CreateString(micm::SolverStateToString(result.state_).c_str()); @@ -253,38 +353,4 @@ namespace musica } } - std::map MICM::GetSpeciesOrdering(Error *error) - { - try - { - micm::State state = solver_->GetState(); - DeleteError(error); - *error = NoError(); - return state.variable_map_; - } - catch (const std::system_error &e) - { - DeleteError(error); - *error = ToError(e); - return std::map(); - } - } - - std::map MICM::GetUserDefinedReactionRatesOrdering(Error *error) - { - try - { - micm::State state = solver_->GetState(); - DeleteError(error); - *error = NoError(); - return state.custom_rate_parameter_map_; - } - catch (const std::system_error &e) - { - DeleteError(error); - *error = ToError(e); - return std::map(); - } - } - } // namespace musica diff --git a/src/test/unit/micm/micm_c_api.cpp b/src/test/unit/micm/micm_c_api.cpp index 3821ceb8..638b7742 100644 --- a/src/test/unit/micm/micm_c_api.cpp +++ b/src/test/unit/micm/micm_c_api.cpp @@ -15,12 +15,13 @@ class MicmCApiTest : public ::testing::Test protected: MICM* micm; const char* config_path = "configs/chapman"; + int num_grid_cells = 1; void SetUp() override { micm = nullptr; Error error; - micm = CreateMicm(config_path, &error); + micm = CreateMicm(config_path, MICMSolver::Rosenbrock, num_grid_cells, &error); ASSERT_TRUE(IsSuccess(error)); DeleteError(&error); @@ -38,13 +39,26 @@ class MicmCApiTest : public ::testing::Test // Test case for bad configuration file path TEST_F(MicmCApiTest, BadConfigurationFilePath) { + int num_grid_cells = 1; Error error = NoError(); - auto micm_bad_config = CreateMicm("bad config path", &error); + auto micm_bad_config = CreateMicm("bad config path", MICMSolver::Rosenbrock, num_grid_cells, &error); ASSERT_EQ(micm_bad_config, nullptr); ASSERT_TRUE(IsError(error, MICM_ERROR_CATEGORY_CONFIGURATION, MICM_CONFIGURATION_ERROR_CODE_INVALID_FILE_PATH)); DeleteError(&error); } +// Test case for bad input for solver type +TEST_F(MicmCApiTest, BadSolverType) +{ + short solver_type = 999; + int num_grid_cells = 1; + Error error = NoError(); + auto micm_bad_solver_type = CreateMicm("configs/chapman", static_cast(solver_type), num_grid_cells, &error); + ASSERT_EQ(micm_bad_solver_type, nullptr); + ASSERT_TRUE(IsError(error, MUSICA_ERROR_CATEGORY, MUSICA_ERROR_CODE_SOLVER_TYPE_NOT_FOUND)); + DeleteError(&error); +} + // Test case for missing species property TEST_F(MicmCApiTest, MissingSpeciesProperty) { @@ -179,8 +193,8 @@ TEST_F(MicmCApiTest, GetUserDefinedReactionRatesOrdering) DeleteMappings(reaction_rates_ordering, array_size); } -// Test case for solving the MICM instance -TEST_F(MicmCApiTest, SolveMicmInstance) +// Test case for solving system using vector-ordered Rosenbrock solver +TEST_F(MicmCApiTest, SolveUsingVectorOrderedRosenbrock) { double time_step = 200.0; double temperature = 272.5; @@ -189,18 +203,19 @@ TEST_F(MicmCApiTest, SolveMicmInstance) double air_density = pressure / (GAS_CONSTANT * temperature); int num_concentrations = 5; double concentrations[] = { 0.75, 0.4, 0.8, 0.01, 0.02 }; + std::size_t num_user_defined_reaction_rates = 3; + double user_defined_reaction_rates[] = { 0.1, 0.2, 0.3 }; String solver_state; SolverResultStats solver_stats; Error error; - auto ordering = micm->GetUserDefinedReactionRatesOrdering(&error); + Mapping* ordering = GetUserDefinedReactionRatesOrdering(micm, &num_user_defined_reaction_rates, &error); ASSERT_TRUE(IsSuccess(error)); - int num_custom_rate_parameters = ordering.size(); - std::vector custom_rate_parameters(num_custom_rate_parameters, 0.0); - for (auto& entry : ordering) + std::vector custom_rate_parameters(num_user_defined_reaction_rates, 0.0); + for (std::size_t i = 0; i < num_user_defined_reaction_rates; i++) { - custom_rate_parameters[entry.second] = 0.0; + custom_rate_parameters[ordering[i].index_] = 0.0; } MicmSolve( @@ -227,7 +242,7 @@ TEST_F(MicmCApiTest, SolveMicmInstance) std::cout << "Solver state: " << solver_state.value_ << std::endl; std::cout << "Function Calls: " << solver_stats.function_calls_ << std::endl; - std::cout << "Jacobian updates:" << solver_stats.jacobian_updates_ << std::endl; + std::cout << "Jacobian updates: " << solver_stats.jacobian_updates_ << std::endl; std::cout << "Number of steps: " << solver_stats.number_of_steps_ << std::endl; std::cout << "Accepted: " << solver_stats.accepted_ << std::endl; std::cout << "Rejected: " << solver_stats.rejected_ << std::endl; @@ -236,10 +251,80 @@ TEST_F(MicmCApiTest, SolveMicmInstance) std::cout << "Singular: " << solver_stats.singular_ << std::endl; std::cout << "Final time: " << solver_stats.final_time_ << std::endl; + DeleteMappings(ordering, num_user_defined_reaction_rates); DeleteString(&solver_state); DeleteError(&error); } +// Test case for solving system using standard-ordered Rosenbrock solver +TEST(RosenbrockStandardOrder, SolveUsingStandardOrderedRosenbrock) +{ + const char* config_path = "configs/chapman"; + int num_grid_cells = 1; + Error error; + MICM* micm = CreateMicm(config_path, MICMSolver::RosenbrockStandardOrder, num_grid_cells, &error); + + double time_step = 200.0; + double temperature = 272.5; + double pressure = 101253.3; + constexpr double GAS_CONSTANT = 8.31446261815324; // J mol-1 K-1 + double air_density = pressure / (GAS_CONSTANT * temperature); + int num_concentrations = 5; + double concentrations[] = { 0.75, 0.4, 0.8, 0.01, 0.02 }; + std::size_t num_user_defined_reaction_rates = 3; + double user_defined_reaction_rates[] = { 0.1, 0.2, 0.3 }; + String solver_state; + SolverResultStats solver_stats; + + Mapping* ordering = GetUserDefinedReactionRatesOrdering(micm, &num_user_defined_reaction_rates, &error); + ASSERT_TRUE(IsSuccess(error)); + + std::vector custom_rate_parameters(num_user_defined_reaction_rates, 0.0); + for (std::size_t i = 0; i < num_user_defined_reaction_rates; i++) + { + custom_rate_parameters[ordering[i].index_] = 0.0; + } + + MicmSolve( + micm, + time_step, + temperature, + pressure, + air_density, + num_concentrations, + concentrations, + custom_rate_parameters.size(), + custom_rate_parameters.data(), + &solver_state, + &solver_stats, + &error); + ASSERT_TRUE(IsSuccess(error)); + + // Add assertions to check the solved concentrations + ASSERT_EQ(concentrations[0], 0.75); + ASSERT_NE(concentrations[1], 0.4); + ASSERT_NE(concentrations[2], 0.8); + ASSERT_NE(concentrations[3], 0.01); + ASSERT_NE(concentrations[4], 0.02); + + std::cout << "Solver state: " << solver_state.value_ << std::endl; + std::cout << "Function Calls: " << solver_stats.function_calls_ << std::endl; + std::cout << "Jacobian updates: " << solver_stats.jacobian_updates_ << std::endl; + std::cout << "Number of steps: " << solver_stats.number_of_steps_ << std::endl; + std::cout << "Accepted: " << solver_stats.accepted_ << std::endl; + std::cout << "Rejected: " << solver_stats.rejected_ << std::endl; + std::cout << "Decompositions: " << solver_stats.decompositions_ << std::endl; + std::cout << "Solves: " << solver_stats.solves_ << std::endl; + std::cout << "Singular: " << solver_stats.singular_ << std::endl; + std::cout << "Final time: " << solver_stats.final_time_ << std::endl; + + DeleteMappings(ordering, num_user_defined_reaction_rates); + DeleteString(&solver_state); + DeleteMicm(micm, &error); + ASSERT_TRUE(IsSuccess(error)); + DeleteError(&error); +} + // Test case for getting species properties TEST_F(MicmCApiTest, GetSpeciesProperty) {