add air density to MICM solve function as input (#160)

NCAR · Jun 24, 2024 · 2d71510 · 2d71510
1 parent e07a447
commit 2d71510
Show file tree

Hide file tree

Showing 8 changed files with 30 additions and 6 deletions.
diff --git a/fortran/micm.F90 b/fortran/micm.F90
@@ -33,14 +33,15 @@ 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
          type(c_ptr), value, intent(in)         :: micm
          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,22 +189,24 @@ 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
       real(c_double),  intent(inout) :: user_defined_reaction_rates(*)
       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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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(),