Add ability to run a PartMC simulation in notebook (#191)

Co-authored-by: Sylwester Arabas <[email protected]>

Author: jcurtis2

.gitmodules

@@ -4,7 +4,7 @@
 	shallow = true
 [submodule "partmc"]
 	path = gitmodules/partmc
-	url = https://github.com/slayoo/partmc
+	url = https://github.com/compdyn/partmc
 	shallow = true
 [submodule "pybind11_json"]
 	path = gitmodules/pybind11_json

.pre-commit-config.yaml

@@ -4,18 +4,18 @@ default_stages: [commit]
 
 repos:
   - repo: https://github.com/psf/black
-    rev: 22.6.0
+    rev: 23.1.0
     hooks:
       - id: black
 
   - repo: https://github.com/timothycrosley/isort
-    rev: 5.10.1
+    rev: 5.12.0
     hooks:
       - id: isort
         args: ["--profile", "black"]
 
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.3.0
+    rev: v4.4.0
     hooks:
       - id: trailing-whitespace
       - id: end-of-file-fixer

CMakeLists.txt

@@ -5,6 +5,11 @@
 ####################################################################################################
 
 cmake_minimum_required(VERSION 3.4...3.18)
+
+if (NOT EXISTS "${CMAKE_SOURCE_DIR}/gitmodules/pybind11/include/pybind11/pybind11.h" )
+  message(FATAL_ERROR "git submodules not initialised.\n Please run `git submodule update --init`")
+endif()
+
 project(_PyPartMC LANGUAGES C CXX Fortran)
 
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)

README.md

@@ -92,10 +92,6 @@ No CMAKE_Fortran_COMPILER could be found.
 ```
 Try installing a Fortran compiler (e.g., `brew reinstall gcc`)
 
-```
-warning: no files found matching 'gitmodules/...
-```
-Since git clone was done without recursive option, try: `git submodule update --init`
 
 ## Notes for developers
 #### How to debug

gitmodules/partmc

@@ -77,7 +77,6 @@ def build_extension(self, ext):  # pylint: disable=too-many-branches
                     pass
 
         else:
-
             # Single config generators are handled "normally"
             single_config = any(x in cmake_generator for x in ("NMake", "Ninja"))
 

setup.py

@@ -42,50 +42,56 @@ subroutine f_aero_state_init(ptr_c, n_part, aero_data_ptr_c) bind(C)
     character(len=10), parameter, dimension(2) :: mode_names = ["init_small","init_large"]
     integer :: n_spec, n_modes, i_mode, i_spec
     integer :: n_part_added, source
+    logical, parameter :: hardcoded_init = .true.
 
     call c_f_pointer(ptr_c, ptr_f)
     call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f)
 
     call aero_state_zero(ptr_f)
-    call aero_state_set_weight(ptr_f, aero_data_ptr_f, AERO_STATE_WEIGHT_FLAT)
-    call aero_state_set_n_part_ideal(ptr_f, n_part)
-
     call fractal_set_spherical(aero_data_ptr_f%fractal)
 
-    n_modes = 2   ! TODO #141
-    n_spec = aero_data_n_spec(aero_data_ptr_f)
-    if (n_spec > 1) then
-    allocate(aero_dist_init%mode(n_modes))
-    do i_mode = 1,n_modes
-       aero_dist_init%mode(i_mode)%name = mode_names(i_mode)
-       aero_dist_init%mode(i_mode)%type = AERO_MODE_TYPE_LOG_NORMAL
-       aero_dist_init%mode(i_mode)%char_radius = diams(i_mode) / 2
-       aero_dist_init%mode(i_mode)%log10_std_dev_radius = & 
-          log10(std(i_mode))
-       aero_dist_init%mode(i_mode)%num_conc = num_conc(i_mode)
-       allocate(aero_dist_init%mode(i_mode)%vol_frac(n_spec))
-       aero_dist_init%mode(i_mode)%vol_frac = 0.0
-       if (i_mode == 1) then
-          i_spec = aero_data_spec_by_name(aero_data_ptr_f, "SO4")
-          aero_dist_init%mode(i_mode)%vol_frac(i_spec) = 1.0
-       else
-          i_spec = aero_data_spec_by_name(aero_data_ptr_f, "OC")
-          aero_dist_init%mode(i_mode)%vol_frac(i_spec) =  .8 
-          i_spec = aero_data_spec_by_name(aero_data_ptr_f, "BC")
-             aero_dist_init%mode(i_mode)%vol_frac(i_spec) =  .2
+    if (hardcoded_init) then
+       n_modes = 2   ! TODO #141
+       n_spec = aero_data_n_spec(aero_data_ptr_f)
+       if (n_spec > 1) then
+          allocate(aero_dist_init%mode(n_modes))
+          do i_mode = 1,n_modes
+             aero_dist_init%mode(i_mode)%name = mode_names(i_mode)
+             aero_dist_init%mode(i_mode)%type = AERO_MODE_TYPE_LOG_NORMAL
+             aero_dist_init%mode(i_mode)%char_radius = diams(i_mode) / 2
+             aero_dist_init%mode(i_mode)%log10_std_dev_radius = &
+                  log10(std(i_mode))
+             aero_dist_init%mode(i_mode)%num_conc = num_conc(i_mode)
+             allocate(aero_dist_init%mode(i_mode)%vol_frac(n_spec))
+             aero_dist_init%mode(i_mode)%vol_frac = 0.0
+             if (i_mode == 1) then
+                i_spec = aero_data_spec_by_name(aero_data_ptr_f, "SO4")
+                aero_dist_init%mode(i_mode)%vol_frac(i_spec) = 1.0
+             else
+                i_spec = aero_data_spec_by_name(aero_data_ptr_f, "OC")
+                aero_dist_init%mode(i_mode)%vol_frac(i_spec) =  .8
+                i_spec = aero_data_spec_by_name(aero_data_ptr_f, "BC")
+                aero_dist_init%mode(i_mode)%vol_frac(i_spec) =  .2
+             end if
+             allocate(aero_dist_init%mode(i_mode)%vol_frac_std(n_spec))
+             aero_dist_init%mode(i_mode)%vol_frac_std = 0.0
+             source = aero_data_source_by_name(aero_data_ptr_f, mode_names(i_mode))
+             aero_dist_init%mode(i_mode)%source = source
+          end do
        end if
-       allocate(aero_dist_init%mode(i_mode)%vol_frac_std(n_spec))
-       aero_dist_init%mode(i_mode)%vol_frac_std = 0.0
-       source = aero_data_source_by_name(aero_data_ptr_f, mode_names(i_mode))
-       aero_dist_init%mode(i_mode)%source = source 
-    end do
-
-    call aero_state_add_aero_dist_sample(ptr_f, aero_data_ptr_f, &
-           aero_dist_init, 1d0, 0d0, &
-           .true., & ! TODO #121 run_part_opt%allow_doubling, &
-           .true., & ! TODO #121 run_part_opt%allow_halving)
-           n_part_added)
     end if
+
+    call aero_state_set_weight(ptr_f, aero_data_ptr_f, AERO_STATE_WEIGHT_NUMMASS_SOURCE)
+    call aero_state_set_n_part_ideal(ptr_f, n_part)
+
+    if (hardcoded_init) then
+       call aero_state_add_aero_dist_sample(ptr_f, aero_data_ptr_f, &
+            aero_dist_init, 1d0, 0d0, &
+            .true., & ! TODO #121 run_part_opt%allow_doubling, &
+            .true., & ! TODO #121 run_part_opt%allow_halving)
+            n_part_added)
+    end if
+
   end subroutine
 
   subroutine f_aero_state_len(ptr_c, len) bind(C)

src/aero_state.F90

@@ -17,6 +17,7 @@ subroutine f_env_state_ctor(ptr_c) bind(C)
         type(c_ptr), intent(out) :: ptr_c
 
         allocate(ptr_f)
+        ptr_f%elapsed_time = 0
         ptr_c = c_loc(ptr_f)
     end subroutine
 
@@ -115,4 +116,26 @@ subroutine f_env_state_get_pressure(ptr_c, pressure) bind(C)
 
     end subroutine
 
+    subroutine f_env_state_get_elapsed_time(ptr_c, elapsed_time) bind(C)
+        type(env_state_t), pointer :: ptr_f => null()
+        type(c_ptr), intent(in) :: ptr_c
+        real(c_double), intent(out) :: elapsed_time 
+
+        call c_f_pointer(ptr_c, ptr_f)
+
+        elapsed_time = ptr_f%elapsed_time
+
+    end subroutine
+
+    subroutine f_env_state_get_start_time(ptr_c, start_time) bind(C)
+        type(env_state_t), pointer :: ptr_f => null()
+        type(c_ptr), intent(in) :: ptr_c
+        real(c_double), intent(out) :: start_time
+
+        call c_f_pointer(ptr_c, ptr_f)
+
+        start_time = ptr_f%start_time
+
+    end subroutine
+
 end module

src/env_state.F90

@@ -19,6 +19,9 @@ extern "C" void f_env_state_set_height(const void *ptr, const double *height) no
 extern "C" void f_env_state_get_height(const void *ptr, double *height) noexcept;
 extern "C" void f_env_state_set_pressure(const void *ptr, const double *pressure) noexcept;
 extern "C" void f_env_state_get_pressure(const void *ptr, double *pressure) noexcept;
+extern "C" void f_env_state_get_elapsed_time(const void *ptr, double *elapsed_time) noexcept;
+extern "C" void f_env_state_get_start_time(const void *ptr, double *start_time) noexcept;
+
 
 struct EnvState {
     PMCResource ptr;
@@ -91,4 +94,24 @@ struct EnvState {
         );
         return pressure;
     }
+
+    static auto get_elapsed_time(const EnvState &self) {
+        double elapsed_time;
+
+        f_env_state_get_elapsed_time(
+            self.ptr.f_arg(),
+            &elapsed_time
+        );
+        return elapsed_time;
+    }
+
+    static auto get_start_time(const EnvState &self) {
+        double start_time;
+
+        f_env_state_get_start_time(
+            self.ptr.f_arg(),
+            &start_time
+        );
+        return start_time;
+    }
 };

src/env_state.hpp

@@ -34,6 +34,8 @@ PYBIND11_MODULE(_PyPartMC, m) {
     )pbdoc";
 
     m.def("run_part", &run_part, "Do a particle-resolved Monte Carlo simulation.");
+    m.def("run_part_timestep", &run_part_timestep, "Do a single time step");
+    m.def("run_part_timeblock", &run_part_timeblock, "Do a time block");
     m.def("condense_equilib_particles", &condense_equilib_particles, R"pbdoc(
       Call condense_equilib_particle() on each particle in the aerosol
       to ensure that every particle has its water content in
@@ -205,6 +207,10 @@ PYBIND11_MODULE(_PyPartMC, m) {
             "returns the current temperature of the environment state")
         .def_property_readonly("rh", EnvState::rh,
             "returns the current relative humidity of the environment state")
+        .def_property_readonly("elapsed_time", EnvState::get_elapsed_time,
+            "returns time since start_time (s).")
+        .def_property_readonly("start_time", EnvState::get_start_time,
+            "returns start time (s since 00:00 UTC on start_day)")
         .def_property("height", &EnvState::get_height, &EnvState::set_height)
         .def_property("pressure", &EnvState::get_pressure, &EnvState::set_pressure)
     ;
@@ -251,6 +257,8 @@ PYBIND11_MODULE(_PyPartMC, m) {
         )
         .def("__str__", Scenario::__str__,
             "returns a string with JSON representation of the object")
+        .def("init_env_state", Scenario::init_env_state,
+            "initializes the EnvState")
     ;
 
     py::class_<GasState>(m,
@@ -386,6 +394,8 @@ PYBIND11_MODULE(_PyPartMC, m) {
         "Scenario",
         "condense_equilib_particles",
         "run_part",
+        "run_part_timeblock",
+        "run_part_timestep",
         "pow2_above",
         "condense_equilib_particle",
         "histogram_1d",