diff --git a/fortran/micm.F90 b/fortran/micm.F90 index 4b59813d..c36ffbd9 100644 --- a/fortran/micm.F90 +++ b/fortran/micm.F90 @@ -33,7 +33,7 @@ subroutine delete_micm_c(micm, error) bind(C, name="DeleteMicm") type(error_t_c), intent(inout) :: error end subroutine delete_micm_c - subroutine micm_solve_c(micm, time_step, temperature, pressure, num_concentrations, concentrations, & + subroutine micm_solve_c(micm, time_step, temperature, pressure, air_density, num_concentrations, concentrations, & num_user_defined_reaction_rates, user_defined_reaction_rates, error) bind(C, name="MicmSolve") use musica_util, only: error_t_c import c_ptr, c_double, c_int @@ -41,6 +41,7 @@ subroutine micm_solve_c(micm, time_step, temperature, pressure, num_concentratio real(kind=c_double), value, intent(in) :: time_step real(kind=c_double), value, intent(in) :: temperature real(kind=c_double), value, intent(in) :: pressure + real(kind=c_double), value, intent(in) :: air_density integer(kind=c_int), value, intent(in) :: num_concentrations real(kind=c_double), intent(inout) :: concentrations(num_concentrations) integer(kind=c_int), value, intent(in) :: num_user_defined_reaction_rates @@ -188,13 +189,14 @@ function constructor(config_path, error) result( this ) end function constructor - subroutine solve(this, time_step, temperature, pressure, num_concentrations, concentrations, & + subroutine solve(this, time_step, temperature, pressure, air_density, num_concentrations, concentrations, & num_user_defined_reaction_rates, user_defined_reaction_rates, error) use musica_util, only: error_t_c, error_t class(micm_t) :: this real(c_double), intent(in) :: time_step real(c_double), intent(in) :: temperature real(c_double), intent(in) :: pressure + real(c_double), intent(in) :: air_density integer(c_int), intent(in) :: num_concentrations real(c_double), intent(inout) :: concentrations(*) integer(c_int), intent(in) :: num_user_defined_reaction_rates @@ -202,8 +204,9 @@ subroutine solve(this, time_step, temperature, pressure, num_concentrations, con type(error_t), intent(out) :: error type(error_t_c) :: error_c - call micm_solve_c(this%ptr, time_step, temperature, pressure, num_concentrations, concentrations, & - num_user_defined_reaction_rates, user_defined_reaction_rates, error_c) + call micm_solve_c(this%ptr, time_step, temperature, pressure, air_density, num_concentrations, & + concentrations, num_user_defined_reaction_rates, user_defined_reaction_rates, & + error_c) error = error_t(error_c) end subroutine solve diff --git a/fortran/test/fetch_content_integration/test_micm_api.F90 b/fortran/test/fetch_content_integration/test_micm_api.F90 index d24b171d..67423a94 100644 --- a/fortran/test/fetch_content_integration/test_micm_api.F90 +++ b/fortran/test/fetch_content_integration/test_micm_api.F90 @@ -27,6 +27,7 @@ subroutine test_api() real(c_double) :: time_step real(c_double) :: temperature real(c_double) :: pressure + real(c_double) :: air_density integer(c_int) :: num_concentrations, num_user_defined_reaction_rates real(c_double), dimension(5) :: concentrations real(c_double), dimension(3) :: user_defined_reaction_rates @@ -37,10 +38,12 @@ subroutine test_api() integer(c_int) :: int_value logical(c_bool) :: bool_value type(error_t) :: error + real(c_double), parameter :: GAS_CONSTANT = 8.31446261815324_c_double ! J mol-1 K-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" @@ -68,7 +71,7 @@ subroutine test_api() 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, & + call micm%solve(time_step, temperature, pressure, air_density, num_concentrations, concentrations, & num_user_defined_reaction_rates, user_defined_reaction_rates, error) ASSERT( error%is_success() ) diff --git a/include/musica/micm.hpp b/include/musica/micm.hpp index 02266906..5dca9e73 100644 --- a/include/musica/micm.hpp +++ b/include/musica/micm.hpp @@ -36,6 +36,7 @@ namespace musica double time_step, double temperature, double pressure, + double air_density, int num_concentrations, double *concentrations, int num_custom_rate_parameters, @@ -63,7 +64,8 @@ namespace musica /// @brief Solve the system /// @param time_step Time [s] to advance the state by /// @param temperature Temperature [K] - /// @param pressure Pressure + /// @param pressure Pressure [Pa] + /// @param air_density Air density [mol m-3] /// @param num_concentrations The size of the concentrations array /// @param concentrations Array of species' concentrations /// @param num_custom_rate_parameters The size of the custom_rate_parameters array @@ -73,6 +75,7 @@ namespace musica double time_step, double temperature, double pressure, + double air_density, int num_concentrations, double *concentrations, int num_custom_rate_parameters, diff --git a/python/test/analytical.py b/python/test/analytical.py index 4dfae1a5..1bd082c4 100644 --- a/python/test/analytical.py +++ b/python/test/analytical.py @@ -8,6 +8,8 @@ def test_simulation(self): time_step = 200.0 temperature = 272.5 pressure = 101253.3 + GAS_CONSTANT = 8.31446261815324 + air_density = pressure / (GAS_CONSTANT * temperature) solver = musica.create_micm("configs/analytical") @@ -47,6 +49,7 @@ def test_simulation(self): time_step, temperature, pressure, + air_density, concentrations, None) model_concentrations.append(concentrations[:]) diff --git a/python/test/chapman.py b/python/test/chapman.py index 4cfebc91..8a16532f 100644 --- a/python/test/chapman.py +++ b/python/test/chapman.py @@ -7,6 +7,8 @@ def test_micm_solve(self): time_step = 200.0 temperature = 272.5 pressure = 101253.3 + GAS_CONSTANT = 8.31446261815324 + air_density = pressure / (GAS_CONSTANT * temperature) concentrations = [0.75, 0.4, 0.8, 0.01, 0.02] solver = musica.create_micm("configs/chapman") @@ -29,6 +31,7 @@ def test_micm_solve(self): time_step, temperature, pressure, + air_density, concentrations, ordered_rate_constants) diff --git a/python/wrapper.cpp b/python/wrapper.cpp index 07778642..0142115d 100644 --- a/python/wrapper.cpp +++ b/python/wrapper.cpp @@ -33,6 +33,7 @@ PYBIND11_MODULE(musica, m) double time_step, double temperature, double pressure, + double air_density, py::list concentrations, py::object custom_rate_parameters = py::none()) { @@ -59,6 +60,7 @@ PYBIND11_MODULE(musica, m) time_step, temperature, pressure, + air_density, concentrations_cpp.size(), concentrations_cpp.data(), custom_rate_parameters_cpp.size(), diff --git a/src/micm/micm.cpp b/src/micm/micm.cpp index 39472577..84e4971f 100644 --- a/src/micm/micm.cpp +++ b/src/micm/micm.cpp @@ -61,6 +61,7 @@ namespace musica double time_step, double temperature, double pressure, + double air_density, int num_concentrations, double *concentrations, int num_custom_rate_parameters, @@ -72,6 +73,7 @@ namespace musica time_step, temperature, pressure, + air_density, num_concentrations, concentrations, num_custom_rate_parameters, @@ -193,6 +195,7 @@ namespace musica double time_step, double temperature, double pressure, + double air_density, int num_concentrations, double *concentrations, int num_custom_rate_parameters, @@ -207,6 +210,7 @@ namespace musica { state.conditions_[i].temperature_ = temperature; state.conditions_[i].pressure_ = pressure; + state.conditions_[i].air_density_ = air_density; } state.variables_.AsVector().assign(concentrations, concentrations + num_concentrations); diff --git a/src/test/unit/micm/micm_c_api.cpp b/src/test/unit/micm/micm_c_api.cpp index e55c0685..49f07bd1 100644 --- a/src/test/unit/micm/micm_c_api.cpp +++ b/src/test/unit/micm/micm_c_api.cpp @@ -183,6 +183,8 @@ TEST_F(MicmCApiTest, SolveMicmInstance) 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 }; @@ -201,6 +203,7 @@ TEST_F(MicmCApiTest, SolveMicmInstance) time_step, temperature, pressure, + air_density, num_concentrations, concentrations, custom_rate_parameters.size(),