diff --git a/.github/workflows/mac.yml b/.github/workflows/mac.yml index 4c99e9e7..247d035b 100644 --- a/.github/workflows/mac.yml +++ b/.github/workflows/mac.yml @@ -29,6 +29,9 @@ jobs: steps: - uses: actions/checkout@v3 + - uses: actions/setup-python@v5 + - run: pip install -r python/requirements.txt + - name: Run Cmake run: cmake -S . -B build -D CMAKE_BUILD_TYPE=${{ matrix.build_type }} -D MUSICA_ENABLE_PYTHON_LIBRARY=ON diff --git a/.github/workflows/ubuntu.yml b/.github/workflows/ubuntu.yml index 4199bab1..b4ad40e4 100644 --- a/.github/workflows/ubuntu.yml +++ b/.github/workflows/ubuntu.yml @@ -29,6 +29,9 @@ jobs: steps: - uses: actions/checkout@v3 + - uses: actions/setup-python@v5 + - run: pip install -r python/requirements.txt + - name: Run Cmake run: cmake -S . -B build -D CMAKE_BUILD_TYPE=${{ matrix.build_type }} -D MUSICA_ENABLE_PYTHON_LIBRARY=ON diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index ee53d20f..f29de4a3 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -25,6 +25,9 @@ jobs: with: platform: ${{ matrix.architecture }} version: 12.2.0 # https://github.com/egor-tensin/setup-mingw/issues/14 + + - uses: actions/setup-python@v5 + - run: pip install -r python/requirements.txt - name: Run Cmake run: cmake -S . -B build -G "MinGW Makefiles" -D MUSICA_ENABLE_PYTHON_LIBRARY=ON @@ -48,11 +51,12 @@ jobs: - uses: actions/checkout@v3 - name: Set up Python - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 if: ${{ matrix.architecture == 'Win32' }} with: python-version: '3.x' architecture: x86 + - run: pip install -r python/requirements.txt - name: Run CMake run: cmake -S . -B build -G "Visual Studio 16 2019" -A ${{ matrix.architecture }} -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} -D MUSICA_ENABLE_PYTHON_LIBRARY=ON @@ -74,12 +78,14 @@ jobs: - uses: actions/checkout@v3 - name: Set up Python - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 if: ${{ matrix.architecture == 'Win32' }} with: python-version: '3.x' architecture: x86 + - run: pip install -r python/requirements.txt + - name: Run CMake run: cmake -S . -B build -G "Visual Studio 17 2022" -A ${{ matrix.architecture }} -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} -D MUSICA_ENABLE_PYTHON_LIBRARY=ON @@ -97,12 +103,19 @@ jobs: steps: - uses: actions/checkout@v3 + - name: Install Clang run: curl -fsSL -o LLVM${{ matrix.version }}.exe https://github.com/llvm/llvm-project/releases/download/llvmorg-${{ matrix.version }}.0.0/LLVM-${{ matrix.version }}.0.0-win64.exe ; 7z x LLVM${{ matrix.version }}.exe -y -o"C:/Program Files/LLVM" + + - uses: actions/setup-python@v5 + - run: pip install -r python/requirements.txt + - name: Run CMake run: cmake -S . -B build -DCMAKE_CXX_COMPILER="C:/Program Files/LLVM/bin/clang++.exe" -DCMAKE_C_COMPILER="C:/Program Files/LLVM/bin/clang.exe" -G"MinGW Makefiles" -D MUSICA_ENABLE_PYTHON_LIBRARY=ON + - name: Build run: cmake --build build --parallel + - name: Test run: cd build ; ctest -j 10 -C Debug --exclude-regex "test-unicode" --output-on-failure @@ -122,6 +135,9 @@ jobs: python-version: '3.x' architecture: x86 + - uses: actions/setup-python@v5 + - run: pip install -r python/requirements.txt + - name: Run CMake run: cmake -S . -B build -G "Visual Studio 16 2019" -A ${{ matrix.architecture }} -T ClangCL -D MUSICA_ENABLE_PYTHON_LIBRARY=ON diff --git a/src/micm/configs/TS1/config.json b/configs/TS1/config.json similarity index 100% rename from src/micm/configs/TS1/config.json rename to configs/TS1/config.json diff --git a/src/micm/configs/TS1/initial_conditions.csv b/configs/TS1/initial_conditions.csv similarity index 100% rename from src/micm/configs/TS1/initial_conditions.csv rename to configs/TS1/initial_conditions.csv diff --git a/src/micm/configs/TS1/reactions.json b/configs/TS1/reactions.json similarity index 100% rename from src/micm/configs/TS1/reactions.json rename to configs/TS1/reactions.json diff --git a/src/micm/configs/TS1/species.json b/configs/TS1/species.json similarity index 100% rename from src/micm/configs/TS1/species.json rename to configs/TS1/species.json diff --git a/configs/analytical/config.json b/configs/analytical/config.json new file mode 100644 index 00000000..03622f15 --- /dev/null +++ b/configs/analytical/config.json @@ -0,0 +1 @@ +{"camp-files": ["species.json", "reactions.json"]} \ No newline at end of file diff --git a/configs/analytical/reactions.json b/configs/analytical/reactions.json new file mode 100644 index 00000000..b44350a3 --- /dev/null +++ b/configs/analytical/reactions.json @@ -0,0 +1 @@ +{"camp-data": [{"type": "MECHANISM", "name": "music box interactive configuration", "reactions": [{"type": "ARRHENIUS", "A": 0.00012, "B": 7, "C" : 75, "D": 50, "E": 0.5, "reactants": {"B": {"qty": 1}}, "products": {"C": {"yield": 1}, "irr__089f1f45-4cd8-4278-83d5-d638e98e4315": {"yield": 1}}}, {"type": "ARRHENIUS", "A": 0.004, "C" : 50, "reactants": {"A": {"qty": 1}}, "products": {"B": {"yield": 1}, "irr__2a109b21-bb24-41ae-8f06-7485fd36f1a7": {"yield": 1}}}]}]} \ No newline at end of file diff --git a/configs/analytical/species.json b/configs/analytical/species.json new file mode 100644 index 00000000..58ede0a7 --- /dev/null +++ b/configs/analytical/species.json @@ -0,0 +1 @@ +{"camp-data": [{"name": "A", "type": "CHEM_SPEC"}, {"name": "B", "type": "CHEM_SPEC"}, {"name": "C", "type": "CHEM_SPEC"}, {"name": "irr__089f1f45-4cd8-4278-83d5-d638e98e4315", "type": "CHEM_SPEC"}, {"name": "irr__2a109b21-bb24-41ae-8f06-7485fd36f1a7", "type": "CHEM_SPEC"}]} \ No newline at end of file diff --git a/src/micm/configs/carbon_bond_5/config.json b/configs/carbon_bond_5/config.json similarity index 100% rename from src/micm/configs/carbon_bond_5/config.json rename to configs/carbon_bond_5/config.json diff --git a/src/micm/configs/carbon_bond_5/initial_conditions.csv b/configs/carbon_bond_5/initial_conditions.csv similarity index 100% rename from src/micm/configs/carbon_bond_5/initial_conditions.csv rename to configs/carbon_bond_5/initial_conditions.csv diff --git a/src/micm/configs/carbon_bond_5/reactions.json b/configs/carbon_bond_5/reactions.json similarity index 100% rename from src/micm/configs/carbon_bond_5/reactions.json rename to configs/carbon_bond_5/reactions.json diff --git a/src/micm/configs/carbon_bond_5/species.json b/configs/carbon_bond_5/species.json similarity index 100% rename from src/micm/configs/carbon_bond_5/species.json rename to configs/carbon_bond_5/species.json diff --git a/configs/chapman/config.json b/configs/chapman/config.json index 04d0ef28..2b5917a0 100644 --- a/configs/chapman/config.json +++ b/configs/chapman/config.json @@ -3,4 +3,4 @@ "species.json", "reactions.json" ] -} +} \ No newline at end of file diff --git a/src/micm/configs/chapman/initial_conditions.csv b/configs/chapman/initial_conditions.csv similarity index 100% rename from src/micm/configs/chapman/initial_conditions.csv rename to configs/chapman/initial_conditions.csv diff --git a/src/micm/configs/robertson/config.json b/configs/robertson/config.json similarity index 100% rename from src/micm/configs/robertson/config.json rename to configs/robertson/config.json diff --git a/src/micm/configs/robertson/initial_conditions.csv b/configs/robertson/initial_conditions.csv similarity index 100% rename from src/micm/configs/robertson/initial_conditions.csv rename to configs/robertson/initial_conditions.csv diff --git a/src/micm/configs/robertson/reactions.json b/configs/robertson/reactions.json similarity index 100% rename from src/micm/configs/robertson/reactions.json rename to configs/robertson/reactions.json diff --git a/src/micm/configs/robertson/species.json b/configs/robertson/species.json similarity index 100% rename from src/micm/configs/robertson/species.json rename to configs/robertson/species.json diff --git a/docker/Dockerfile.fortran-intel b/docker/Dockerfile.fortran-intel index 8e90f135..e3bd63d7 100644 --- a/docker/Dockerfile.fortran-intel +++ b/docker/Dockerfile.fortran-intel @@ -33,7 +33,7 @@ RUN cd musica \ -B build \ -D MUSICA_ENABLE_TESTS=ON \ -D MUSICA_ENABLE_TUVX=OFF \ - -D CMAKE_BUILD_TYPE=Debug \ + -D CMAKE_BUILD_TYPE=Release \ && cd build \ && make install @@ -45,7 +45,7 @@ ENV FC=ifx RUN cd musica/fortran/test/fetch_content_integration \ && mkdir build && cd build \ && cmake .. \ - -D CMAKE_BUILD_TYPE=Debug \ + -D CMAKE_BUILD_TYPE=Release \ && make RUN cp -r /musica/configs/chapman musica/fortran/test/fetch_content_integration/build diff --git a/docker/Dockerfile.fortran-nvhpc b/docker/Dockerfile.fortran-nvhpc index 9c21239a..b86296b8 100644 --- a/docker/Dockerfile.fortran-nvhpc +++ b/docker/Dockerfile.fortran-nvhpc @@ -38,7 +38,7 @@ RUN cd musica \ -B build \ -D MUSICA_ENABLE_TESTS=ON \ -D MUSICA_ENABLE_TUVX=OFF \ - -D CMAKE_BUILD_TYPE=Debug \ + -D CMAKE_BUILD_TYPE=Release \ && cd build \ && make install @@ -50,7 +50,7 @@ ENV FC=nvfortran RUN cd musica/fortran/test/fetch_content_integration \ && mkdir build && cd build \ && cmake .. \ - -D CMAKE_BUILD_TYPE=Debug \ + -D CMAKE_BUILD_TYPE=Release \ && make RUN cp -r /musica/configs/chapman musica/fortran/test/fetch_content_integration/build diff --git a/fortran/micm_core.F90 b/fortran/micm_core.F90 index 80d3cc18..35760999 100644 --- a/fortran/micm_core.F90 +++ b/fortran/micm_core.F90 @@ -1,11 +1,16 @@ module micm_core - use iso_c_binding, only: c_ptr, c_char, c_int, c_double, c_null_char + use iso_c_binding, only: c_ptr, c_char, c_int, c_double, c_null_char, c_size_t, c_f_pointer implicit none public :: micm_t private + type, bind(c) :: Mapping + character(kind=c_char, len=1) :: name(256) + integer(c_size_t) :: index + end type Mapping + interface function create_micm_c(config_path, error_code) bind(C, name="create_micm") import c_ptr, c_int, c_char @@ -28,11 +33,26 @@ subroutine micm_solve_c(micm, time_step, temperature, pressure, num_concentratio integer(kind=c_int), value, intent(in) :: num_concentrations real(kind=c_double), intent(inout) :: concentrations(num_concentrations) end subroutine micm_solve_c + + function get_species_ordering(micm, array_size) bind(c, name="get_species_ordering") + import :: c_ptr, c_size_t + type(c_ptr), value :: micm + integer(kind=c_size_t), intent(out) :: array_size + type(c_ptr) :: get_species_ordering + end function get_species_ordering + + type(c_ptr) function get_user_defined_reaction_rates_ordering(micm, array_size) & + bind(c, name="get_user_defined_reaction_rates_ordering") + import :: c_ptr, c_size_t + type(c_ptr), value :: micm + integer(kind=c_size_t), intent(out) :: array_size + end function get_user_defined_reaction_rates_ordering end interface type :: micm_t - private - type(c_ptr) :: ptr + type(Mapping), pointer :: species_ordering(:), reaction_rates_ordering(:) + integer(kind=c_size_t) :: species_ordering_length, reaction_rates_ordering_length + type(c_ptr), private :: ptr contains ! Solve the chemical system procedure :: solve @@ -52,9 +72,11 @@ function constructor(config_path, errcode) result( this ) integer, intent(out) :: errcode character(len=1, kind=c_char) :: c_config_path(len_trim(config_path)+1) integer :: n, i + type(c_ptr) :: mappings_ptr allocate( this ) + print *, "Config Path:", config_path n = len_trim(config_path) do i = 1, n c_config_path(i) = config_path(i:i) @@ -66,6 +88,15 @@ function constructor(config_path, errcode) result( this ) if (errcode /= 0) then return end if + + mappings_ptr = get_species_ordering(this%ptr, this%species_ordering_length) + call c_f_pointer(mappings_ptr, this%species_ordering, [this%species_ordering_length]) + + ! this%reaction_rates_ordering = get_user_defined_reaction_rates_ordering(this%ptr, this%reaction_rates_ordering_length) + + ! do i = 1, this%reaction_rates_ordering_length + ! print *, "Reaction Rate Name:", this%reaction_rates_ordering(i)%name, ", Index:", this%reaction_rates_ordering(i)%index + ! end do end function constructor subroutine solve(this, time_step, temperature, pressure, num_concentrations, concentrations) diff --git a/fortran/test/fetch_content_integration/CMakeLists.txt b/fortran/test/fetch_content_integration/CMakeLists.txt index e5e00d46..10ffe012 100644 --- a/fortran/test/fetch_content_integration/CMakeLists.txt +++ b/fortran/test/fetch_content_integration/CMakeLists.txt @@ -11,7 +11,7 @@ include(FetchContent) FetchContent_Declare(musica-fortran GIT_REPOSITORY https://github.com/NCAR/musica.git - GIT_TAG main + GIT_TAG 79-allow-updates-to-custom-rate-constant-parameters-and-mapping-by-name ) set(MUSICA_BUILD_C_CXX_INTERFACE OFF) diff --git a/fortran/test/unit/CMakeLists.txt b/fortran/test/unit/CMakeLists.txt index f67af53f..4e3062c2 100644 --- a/fortran/test/unit/CMakeLists.txt +++ b/fortran/test/unit/CMakeLists.txt @@ -1,5 +1,7 @@ include(test_util) +create_standard_test_fortran(NAME musica_fortran_interface SOURCES test_musica_fortran_interface.F90) + if (MUSICA_ENABLE_TUVX) create_standard_test_fortran(NAME connect_to_tuvx SOURCES tuvx.F90) if (MUSICA_ENABLE_OPENMP) diff --git a/fortran/test/unit/test_musica_fortran_interface.F90 b/fortran/test/unit/test_musica_fortran_interface.F90 new file mode 100644 index 00000000..f7d6ca2e --- /dev/null +++ b/fortran/test/unit/test_musica_fortran_interface.F90 @@ -0,0 +1,40 @@ +program test_micm_fort_api + use iso_c_binding + use micm_core, only: micm_t + + implicit none + + type(micm_t), pointer :: micm + real(c_double) :: time_step + real(c_double) :: temperature + real(c_double) :: pressure + integer(c_int) :: num_concentrations + real(c_double), dimension(5) :: concentrations + integer :: errcode + character(len=7) :: config_path + + time_step = 200 + temperature = 272.5 + pressure = 101253.3 + num_concentrations = 5 + concentrations = (/ 0.75, 0.4, 0.8, 0.01, 0.02 /) + config_path = "configs/chapman" + + write(*,*) "[test micm fort api] Creating MICM solver..." + micm => micm_t(config_path, errcode) + + if (errcode /= 0) then + write(*,*) "[test micm fort api] Failed in creating solver." + stop 3 + endif + + write(*,*) "[test micm fort api] Initial concentrations", concentrations + + write(*,*) "[test micm fort api] Solving starts..." + call micm%solve(time_step, temperature, pressure, num_concentrations, concentrations) + + write(*,*) "[test micm fort api] After solving, concentrations", concentrations + + write(*,*) "[test micm fort api] Finished." + +end program diff --git a/include/musica/micm.hpp b/include/musica/micm.hpp index 4a680363..e4418df9 100644 --- a/include/musica/micm.hpp +++ b/include/musica/micm.hpp @@ -2,7 +2,7 @@ * This file contains the defintion of the MICM class, which represents a multi-component * reactive transport model. It also includes functions for creating and deleting MICM instances with c binding * Copyright (C) 2023-2024 National Center for Atmospheric Research, - * + * * SPDX-License-Identifier: Apache-2.0* creating solvers, and solving the model. */ #pragma once @@ -15,13 +15,22 @@ class MICM; +struct Mapping +{ + char *name; + size_t index; +}; + #ifdef __cplusplus -extern "C" { +extern "C" +{ #endif -MICM* create_micm(const char* config_path, int* error_code); -void delete_micm(const MICM* micm); -void micm_solve(MICM* micm, double time_step, double temperature, double pressure, int num_concentrations, double* concentrations); + MICM *create_micm(const char *config_path, int *error_code); + void delete_micm(const MICM *micm); + void micm_solve(MICM *micm, double time_step, double temperature, double pressure, int num_concentrations, double *concentrations, int num_custom_rate_parameters, double *custom_rate_parameters); + Mapping *get_species_ordering(MICM *micm, size_t *array_size); + Mapping *get_user_defined_reaction_rates_ordering(MICM *micm, size_t *array_size); #ifdef __cplusplus } @@ -39,8 +48,8 @@ class MICM /// @brief Create a solver by reading and parsing configuration file /// @param config_path Path to configuration file or directory containing configuration file /// @return 0 on success, 1 on failure in parsing configuration file - int create_solver(const std::string& config_path); - + int create_solver(const std::string &config_path); + /// @brief Solve the system /// @param time_step Time [s] to advance the state by /// @para/ @briefm temperature Temperature [K] @@ -48,7 +57,10 @@ class MICM /// @param config_path Path to configuration file or directory containing configuration file /// @return 0 on success, 1 on failure in parsing species' concentrations /// @param concentrations Species's concentrations - void solve(double time_step, double temperature, double pressure, int num_concentrations, double* concentrations); + void solve(double time_step, double temperature, double pressure, int num_concentrations, double *concentrations, int num_custom_rate_parameters, double *custom_rate_parameters); + + std::map get_species_ordering(); + std::map get_user_defined_reaction_rates_ordering(); static constexpr size_t NUM_GRID_CELLS = 1; diff --git a/lib/micm b/lib/micm index 4a9149db..9993cd95 160000 --- a/lib/micm +++ b/lib/micm @@ -1 +1 @@ -Subproject commit 4a9149db80f4f5adfb4cb751b70b0f7818f23616 +Subproject commit 9993cd953aa33e96d677c6ba532e6e660bc0a43b diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 4b14208c..0b457c15 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -12,7 +12,6 @@ target_link_libraries(musica_python PUBLIC musica::musica) if(WIN32) - # makefiles on windows don't need the config directory if (${CMAKE_GENERATOR} MATCHES "MinGW Makefiles") set(PYTHON_MODULE_PATH "${CMAKE_CURRENT_BINARY_DIR}") @@ -24,7 +23,5 @@ else() set(PYTHON_MODULE_PATH "${CMAKE_CURRENT_BINARY_DIR}") endif() -add_test(NAME musica_python_test - WORKING_DIRECTORY ${CMAKE_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env PYTHONPATH=${PYTHON_MODULE_PATH} - ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/micm_test.py) \ No newline at end of file + +add_subdirectory(test) \ No newline at end of file diff --git a/python/micm_test.py b/python/micm_test.py deleted file mode 100644 index de926c78..00000000 --- a/python/micm_test.py +++ /dev/null @@ -1,23 +0,0 @@ -import sys - -import musica - -#tests python micm package -time_step = 200.0 -temperature = 272.5 -pressure = 101253.3 -num_concentrations = 5 -concentrations = [0.75, 0.4, 0.8, 0.01, 0.02] - -print(concentrations) - -solver = musica.create_micm("configs/chapman") -musica.micm_solve(solver, time_step, temperature, pressure, concentrations) - -print(concentrations) - -assert concentrations[0] == 0.75 -assert concentrations[1] != 0.4 -assert concentrations[2] != 0.8 -assert concentrations[3] != 0.01 -assert concentrations[4] != 0.02 diff --git a/python/requirements.txt b/python/requirements.txt new file mode 100644 index 00000000..296d6545 --- /dev/null +++ b/python/requirements.txt @@ -0,0 +1 @@ +numpy \ No newline at end of file diff --git a/python/test/CMakeLists.txt b/python/test/CMakeLists.txt new file mode 100644 index 00000000..04fd7ac9 --- /dev/null +++ b/python/test/CMakeLists.txt @@ -0,0 +1,9 @@ +add_test(NAME musica_python_chapman_test + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + COMMAND ${CMAKE_COMMAND} -E env PYTHONPATH=${PYTHON_MODULE_PATH} + ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/chapman.py) + +add_test(NAME musica_python_analytical_test + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + COMMAND ${CMAKE_COMMAND} -E env PYTHONPATH=${PYTHON_MODULE_PATH} + ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/analytical.py) \ No newline at end of file diff --git a/python/test/analytical.py b/python/test/analytical.py new file mode 100644 index 00000000..685e2519 --- /dev/null +++ b/python/test/analytical.py @@ -0,0 +1,60 @@ +import unittest +import numpy as np +import musica + +class TestAnalyticalSimulation(unittest.TestCase): + def test_simulation(self): + time_step = 200.0 + temperature = 272.5 + pressure = 101253.3 + + solver = musica.create_micm("configs/analytical") + + rates = musica.user_defined_reaction_rates(solver) + ordering = musica.species_ordering(solver) + + self.assertEqual(rates, {}) + self.assertEqual(ordering['A'], 0) + self.assertEqual(ordering['B'], 1) + self.assertEqual(ordering['C'], 2) + + # Initalizes concentrations + initial_concentrations = [1, 0, 0] + + model_concentrations = [initial_concentrations[:]] + analytical_concentrations = [initial_concentrations[:]] + + k1 = 4.0e-3 * np.exp(50 / temperature) + k2 = 1.2e-4 * np.exp(75 / temperature) * (temperature / 50)**7* (1.0 + 0.5 * pressure) + + time_step = 1 + sim_length = 100 + + curr_time = time_step + + idx_A = 0 + idx_B = 1 + idx_C = 2 + + concentrations = initial_concentrations + + # Gets analytical concentrations + while curr_time <= sim_length: + musica.micm_solve(solver, time_step, temperature, pressure, concentrations, None) + model_concentrations.append(concentrations[:]) + + initial_A = analytical_concentrations[0][idx_A] + A_conc = initial_A * np.exp(-(k1)* curr_time) + B_conc = initial_A * (k1 / (k2 - k1)) * (np.exp(-k1 * curr_time) - np.exp(-k2 * curr_time)) + C_conc = initial_A * (1.0 + (k1 * np.exp(-k2 * curr_time) - k2 * np.exp(-k1 * curr_time)) / (k2 - k1)) + + analytical_concentrations.append([A_conc, B_conc, C_conc]) + curr_time += time_step + + for i in range(len(model_concentrations)): + self.assertAlmostEqual(model_concentrations[i][idx_A], analytical_concentrations[i][idx_A], places=8) + self.assertAlmostEqual(model_concentrations[i][idx_B], analytical_concentrations[i][idx_B], places=8) + self.assertAlmostEqual(model_concentrations[i][idx_C], analytical_concentrations[i][idx_C], places=8) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/python/test/chapman.py b/python/test/chapman.py new file mode 100644 index 00000000..87b46ec8 --- /dev/null +++ b/python/test/chapman.py @@ -0,0 +1,38 @@ +import unittest +import musica + +class TestChapman(unittest.TestCase): + def test_micm_solve(self): + time_step = 200.0 + temperature = 272.5 + pressure = 101253.3 + concentrations = [0.75, 0.4, 0.8, 0.01, 0.02] + + solver = musica.create_micm("configs/chapman") + rate_constant_ordering = musica.user_defined_reaction_rates(solver) + ordering = musica.species_ordering(solver) + + rate_constants = { + "PHOTO.R1" : 2.42e-17, + "PHOTO.R3" : 1.15e-5, + "PHOTO.R5" : 6.61e-9 + } + + ordered_rate_constants = len(rate_constants.keys()) * [0.0] + + for key, value in rate_constants.items(): + ordered_rate_constants[rate_constant_ordering[key]] = value + + musica.micm_solve(solver, time_step, temperature, pressure, concentrations, ordered_rate_constants) + + self.assertEqual(ordering, {'M': 0, 'O': 2, 'O1D': 1, 'O2': 3, 'O3': 4}) + self.assertEqual(rate_constant_ordering, {'PHOTO.R1': 0, 'PHOTO.R3': 1, 'PHOTO.R5': 2}) + + self.assertEqual(concentrations[0], 0.75) + self.assertNotEqual(concentrations[1], 0.4) + self.assertNotEqual(concentrations[2], 0.8) + self.assertNotEqual(concentrations[3], 0.01) + self.assertNotEqual(concentrations[4], 0.02) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/python/wrapper.cpp b/python/wrapper.cpp index e3390671..b1cec1a8 100644 --- a/python/wrapper.cpp +++ b/python/wrapper.cpp @@ -1,38 +1,59 @@ #include -#include +#include #include namespace py = pybind11; -//Wraps micm.cpp -PYBIND11_MODULE(musica, m) { +// Wraps micm.cpp +PYBIND11_MODULE(musica, m) +{ py::class_(m, "MICM") .def(py::init<>()) .def("create_solver", &MICM::create_solver) .def("solve", &MICM::solve) - .def("__del__", [](MICM &micm) { - std::cout << "MICM destructor called" << std::endl; - }); + .def("__del__", [](MICM &micm) {}); - m.def("create_micm", [](const char* config_path) { + m.def("create_micm", [](const char *config_path) + { int error_code; MICM* micm = create_micm(config_path, &error_code); - return micm; - }); + return micm; }); m.def("delete_micm", &delete_micm); - m.def("micm_solve", [](MICM* micm, double time_step, double temperature, double pressure, py::list concentrations) { + m.def( + "micm_solve", [](MICM *micm, double time_step, double temperature, double pressure, py::list concentrations, py::object custom_rate_parameters = py::none()) + { std::vector concentrations_cpp; for (auto item : concentrations) { concentrations_cpp.push_back(item.cast()); } - micm_solve(micm, time_step, temperature, pressure, concentrations_cpp.size(), concentrations_cpp.data()); + + std::vector custom_rate_parameters_cpp; + if (!custom_rate_parameters.is_none()) { + py::list parameters = custom_rate_parameters.cast(); + for (auto item : parameters) { + custom_rate_parameters_cpp.push_back(item.cast()); + } + } + + micm_solve(micm, time_step, temperature, pressure, + concentrations_cpp.size(), concentrations_cpp.data(), + custom_rate_parameters_cpp.size(), custom_rate_parameters_cpp.data()); // Update the concentrations list after solving for (size_t i = 0; i < concentrations_cpp.size(); ++i) { concentrations[i] = concentrations_cpp[i]; - } - }); + } }, + "Solve the system"); + + m.def( + "species_ordering", [](MICM *micm) + { return micm->get_species_ordering(); }, + "Return map of get_species_ordering rates"); + m.def( + "user_defined_reaction_rates", [](MICM *micm) + { return micm->get_user_defined_reaction_rates_ordering(); }, + "Return map of reaction rates"); } \ No newline at end of file diff --git a/src/micm/__pycache__/micm.cpython-39.pyc b/src/micm/__pycache__/micm.cpython-39.pyc deleted file mode 100644 index d191ba5c..00000000 Binary files a/src/micm/__pycache__/micm.cpython-39.pyc and /dev/null differ diff --git a/src/micm/configs/chapman/config.json b/src/micm/configs/chapman/config.json deleted file mode 100644 index 04d0ef28..00000000 --- a/src/micm/configs/chapman/config.json +++ /dev/null @@ -1,6 +0,0 @@ -{ - "camp-files": [ - "species.json", - "reactions.json" - ] -} diff --git a/src/micm/configs/chapman/reactions.json b/src/micm/configs/chapman/reactions.json deleted file mode 100644 index 679593f9..00000000 --- a/src/micm/configs/chapman/reactions.json +++ /dev/null @@ -1,97 +0,0 @@ -{ - "camp-data": [ - { - "name": "Chapman", - "type": "MECHANISM", - "reactions": [ - { - "type": "PHOTOLYSIS", - "reactants": { - "O2": {} - }, - "products": { - "O": { - "yield": 2.0 - } - }, - "MUSICA name": "R1" - }, - { - "type": "ARRHENIUS", - "A": 8.018e-17, - "reactants": { - "O": {}, - "O2": {} - }, - "products": { - "O3": {} - }, - "MUSICA name": "R2" - }, - { - "type": "PHOTOLYSIS", - "reactants": { - "O3": {} - }, - "products": { - "O": {}, - "O2": {} - }, - "MUSICA name": "R3" - }, - { - "type": "ARRHENIUS", - "A": 1.576e-15, - "reactants": { - "O": {}, - "O3": {} - }, - "products": { - "O2": { - "yield": 2.0 - } - }, - "MUSICA name": "R4" - }, - { - "type": "PHOTOLYSIS", - "reactants": { - "O3": {} - }, - "products": { - "O1D": {}, - "O2": {} - }, - "MUSICA name": "R5" - }, - { - "type": "ARRHENIUS", - "A": 7.11e-11, - "reactants": { - "O1D": {}, - "M": {} - }, - "products": { - "O": {}, - "M": {} - }, - "MUSICA name": "R6" - }, - { - "type": "ARRHENIUS", - "A": 1.2e-10, - "reactants": { - "O1D": {}, - "O3": {} - }, - "products": { - "O2": { - "yield": 2.0 - } - }, - "MUSICA name": "R7" - } - ] - } - ] -} \ No newline at end of file diff --git a/src/micm/configs/chapman/species.json b/src/micm/configs/chapman/species.json deleted file mode 100644 index c962907f..00000000 --- a/src/micm/configs/chapman/species.json +++ /dev/null @@ -1,29 +0,0 @@ -{ - "camp-data": [ - { - "name": "M", - "type": "CHEM_SPEC", - "tracer type": "CONSTANT" - }, - { - "name": "O2", - "type": "CHEM_SPEC", - "tracer type": "CONSTANT" - }, - { - "name": "O", - "type": "CHEM_SPEC", - "absolute tolerance": 1e-12 - }, - { - "name": "O1D", - "type": "CHEM_SPEC", - "absolute tolerance": 1e-12 - }, - { - "name": "O3", - "type": "CHEM_SPEC", - "absolute tolerance": 1e-12 - } - ] -} \ No newline at end of file diff --git a/src/micm/micm.cpp b/src/micm/micm.cpp index f8d6fa19..8dded4cc 100644 --- a/src/micm/micm.cpp +++ b/src/micm/micm.cpp @@ -2,7 +2,7 @@ * This file contains the implementation of the MICM class, which represents a multi-component * reactive transport model. It also includes functions for creating and deleting MICM instances, * Copyright (C) 2023-2024 National Center for Atmospheric Research, - * + * * SPDX-License-Identifier: Apache-2.0* creating solvers, and solving the model. */ #include @@ -12,27 +12,71 @@ #include #include -MICM* create_micm(const char* config_path, int* error_code) +MICM *create_micm(const char *config_path, int *error_code) { - try { - MICM* micm = new MICM(); + try + { + MICM *micm = new MICM(); *error_code = micm->create_solver(std::string(config_path)); return micm; } - catch (const std::bad_alloc& e) { + catch (const std::bad_alloc &e) + { *error_code = 1; return nullptr; } } -void delete_micm(const MICM* micm) +void delete_micm(const MICM *micm) { delete micm; } -void micm_solve(MICM* micm, double time_step, double temperature, double pressure, int num_concentrations, double* concentrations) +void micm_solve(MICM *micm, double time_step, double temperature, double pressure, int num_concentrations, double *concentrations, int num_custom_rate_parameters, double *custom_rate_parameters) +{ + micm->solve(time_step, temperature, pressure, num_concentrations, concentrations, num_custom_rate_parameters, custom_rate_parameters); +} + +Mapping *get_species_ordering(MICM *micm, size_t *array_size) +{ + auto map = micm->get_species_ordering(); + Mapping *reactionRates = new Mapping[map.size()]; + + // Copy data from the map to the array of structs + size_t i = 0; + for (const auto &entry : map) + { + reactionRates[i].name = new char[entry.first.size() + 1]; // +1 for null terminator + std::strcpy(reactionRates[i].name, entry.first.c_str()); + reactionRates[i].index = entry.second; + ++i; + } + + // Set the size of the array + *array_size = map.size(); + + return reactionRates; +} + +Mapping *get_user_defined_reaction_rates_ordering(MICM *micm, size_t *array_size) { - micm->solve(time_step, temperature, pressure, num_concentrations, concentrations); + auto map = micm->get_user_defined_reaction_rates_ordering(); + Mapping *reactionRates = new Mapping[map.size()]; + + // Copy data from the map to the array of structs + size_t i = 0; + for (const auto &entry : map) + { + reactionRates[i].name = new char[entry.first.size() + 1]; // +1 for null terminator + std::strcpy(reactionRates[i].name, entry.first.c_str()); + reactionRates[i].index = entry.second; + ++i; + } + + // Set the size of the array + *array_size = map.size(); + + return reactionRates; } MICM::MICM() : solver_(nullptr) {} @@ -44,27 +88,35 @@ MICM::~MICM() int MICM::create_solver(const std::string &config_path) { int parsing_status = 0; // 0 on success, 1 on failure + try { + micm::SolverConfig<> solver_config; + micm::ConfigParseStatus status = solver_config.ReadAndParse(std::filesystem::path(config_path)); - micm::SolverConfig<> solver_config; - micm::ConfigParseStatus status = solver_config.ReadAndParse(std::filesystem::path(config_path)); - - if (status == micm::ConfigParseStatus::Success) - { - micm::SolverParameters solver_params = solver_config.GetSolverParams(); - auto params = micm::RosenbrockSolverParameters::three_stage_rosenbrock_parameters(NUM_GRID_CELLS); - solver_ = std::make_unique>(solver_params.system_, - solver_params.processes_, - params); + if (status == micm::ConfigParseStatus::Success) + { + micm::SolverParameters solver_params = solver_config.GetSolverParams(); + auto params = micm::RosenbrockSolverParameters::three_stage_rosenbrock_parameters(NUM_GRID_CELLS); + params.ignore_unused_species_ = true; + solver_ = std::make_unique>(solver_params.system_, + solver_params.processes_, + params); + } + else + { + parsing_status = 1; + } } - else + catch(std::exception &e) { + std::cerr << "Error: " << e.what() << std::endl; parsing_status = 1; } + return parsing_status; } -void MICM::solve(double time_step, double temperature, double pressure, int num_concentrations, double *concentrations) +void MICM::solve(double time_step, double temperature, double pressure, int num_concentrations, double *concentrations, int num_custom_rate_parameters, double *custom_rate_parameters) { micm::State state = solver_->GetState(); @@ -76,6 +128,8 @@ void MICM::solve(double time_step, double temperature, double pressure, int num_ state.variables_.AsVector().assign(concentrations, concentrations + num_concentrations); + state.custom_rate_parameters_.AsVector().assign(custom_rate_parameters, custom_rate_parameters + num_custom_rate_parameters); + auto result = solver_->Solve(time_step, state); for (int i = 0; i < result.result_.AsVector().size(); i++) @@ -83,3 +137,15 @@ void MICM::solve(double time_step, double temperature, double pressure, int num_ concentrations[i] = result.result_.AsVector()[i]; } } + +std::map MICM::get_species_ordering() +{ + micm::State state = solver_->GetState(); + return state.variable_map_; +} + +std::map MICM::get_user_defined_reaction_rates_ordering() +{ + micm::State state = solver_->GetState(); + return state.custom_rate_parameter_map_; +} diff --git a/src/test/connections/micm_c_api.cpp b/src/test/connections/micm_c_api.cpp index c7ee270c..a58fce23 100644 --- a/src/test/connections/micm_c_api.cpp +++ b/src/test/connections/micm_c_api.cpp @@ -33,7 +33,14 @@ TEST_F(MicmCApiTest, SolveMicmInstance) { int num_concentrations = 5; double concentrations[] = {0.75, 0.4, 0.8, 0.01, 0.02}; - micm_solve(micm, time_step, temperature, pressure, num_concentrations, concentrations); + auto ordering = micm->get_user_defined_reaction_rates_ordering(); + int num_custom_rate_parameters = ordering.size(); + std::vector custom_rate_parameters(num_custom_rate_parameters, 0.0); + for(auto& entry : ordering) { + custom_rate_parameters[entry.second] = 0.0; + } + + micm_solve(micm, time_step, temperature, pressure, num_concentrations, concentrations, custom_rate_parameters.size(), custom_rate_parameters.data()); // Add assertions to check the solved concentrations ASSERT_EQ(concentrations[0], 0.75);