remove hardcoded AeroState init parameters, add AeroState::dist_sample. closes #116 (#197)

Co-authored-by: Jeff Curtis <[email protected]>

Author: slayoo

src/aero_state.F90

@@ -33,65 +33,13 @@ subroutine f_aero_state_init(ptr_c, n_part, aero_data_ptr_c) bind(C)
     type(c_ptr) :: ptr_c, aero_data_ptr_c
     real(c_double), intent(in) :: n_part
 
-    ! local testing ! TODO #141
-    type(aero_dist_t) :: aero_dist_init
-    real(kind=dp), parameter, dimension(2) :: num_conc = &
-         [3.2d9, 2.9d9]  ! TODO #141
-    real(kind=dp), parameter, dimension(2) :: diams = [2d-8, 1.16d-7]  ! TODO #141
-    real(kind=dp), parameter, dimension(2) :: std = [1.45,1.65]  ! TODO #141
-    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 fractal_set_spherical(aero_data_ptr_f%fractal)
-
-    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
-    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)
@@ -238,7 +186,7 @@ subroutine f_aero_state_crit_rel_humids(ptr_c, aero_data_ptr_c, &
 
   ! TODO #130: Add include, exclude, group and groups
   subroutine f_aero_state_mixing_state_metrics(ptr_c, aero_data_ptr_c, & 
-       d_alpha, d_gamma, chi) bind(C) !, include, exclude, group, groups)&
+       d_alpha, d_gamma, chi) bind(C)
 
     type(aero_state_t), pointer :: ptr_f => null()
     type(aero_data_t), pointer :: aero_data_ptr_f => null()
@@ -310,4 +258,26 @@ subroutine f_aero_state_rand_particle(ptr_c, ptr_particle_c) bind(C)
 
   end subroutine
 
+  subroutine f_aero_state_add_aero_dist_sample(ptr_c, ptr_aero_data_c, &
+       ptr_aero_dist_c, sample_prop, create_time, allow_doubling, &
+       allow_halving, n_part_add) bind(C)
+
+    type(c_ptr) :: ptr_c, ptr_aero_data_c, ptr_aero_dist_c
+    type(aero_state_t), pointer :: ptr_f => null()
+    type(aero_data_t), pointer :: ptr_aero_data_f => null()
+    type(aero_dist_t), pointer :: ptr_aero_dist_f => null()
+    real(c_double) :: sample_prop, create_time
+    logical(c_bool) :: allow_doubling, allow_halving
+    integer(c_int) :: n_part_add
+
+    call c_f_pointer(ptr_c, ptr_f)
+    call c_f_pointer(ptr_aero_data_c,ptr_aero_data_f)
+    call c_f_pointer(ptr_aero_dist_c,ptr_aero_dist_f)
+
+    call aero_state_add_aero_dist_sample(ptr_f, ptr_aero_data_f, &
+       ptr_aero_dist_f, sample_prop, create_time, LOGICAL(allow_doubling), &
+       logical(allow_halving), n_part_add)
+
+  end subroutine
+
 end module

src/aero_state.hpp

@@ -8,6 +8,7 @@
 
 #include "pmc_resource.hpp"
 #include "aero_data.hpp"
+#include "aero_dist.hpp"
 #include "aero_particle.hpp"
 #include "env_state.hpp"
 #include "bin_grid.hpp"
@@ -116,6 +117,17 @@ extern "C" void f_aero_state_rand_particle(
     const void *ptr_particle_c
 ) noexcept;
 
+extern "C" void f_aero_state_add_aero_dist_sample(
+    const void *ptr_c, 
+    const void *ptr_aero_data_c, 
+    const void *ptr_aero_dist_c, 
+    const double *sample_prop,
+    const double *create_time,
+    const bool *allow_doubling,
+    const bool *allow_halving,
+    int *n_part_add
+) noexcept;
+
 struct AeroState {
     PMCResource ptr;
     std::shared_ptr<AeroData> aero_data;
@@ -355,4 +367,28 @@ struct AeroState {
 
         return ptr;
     }
+
+   static int dist_sample(
+       const AeroState &self,
+       const AeroDist &aero_dist,
+       const double &sample_prop,
+       const double &create_time,
+       const bool &allow_doubling,
+       const bool &allow_halving
+   ) {
+       int n_part_add = 0;
+
+       f_aero_state_add_aero_dist_sample(
+       self.ptr.f_arg(),
+       self.aero_data->ptr.f_arg(),
+       aero_dist.ptr.f_arg(),
+       &sample_prop,
+       &create_time,
+       &allow_doubling,
+       &allow_halving,
+       &n_part_add
+       );
+       return n_part_add;
+
+   }
 };

src/pypartmc.cpp

@@ -167,6 +167,8 @@ PYBIND11_MODULE(_PyPartMC, m) {
             "returns the particle of a given index")
         .def("rand_particle", AeroState::get_random_particle,
             "returns a random particle from the population")
+        .def("dist_sample", AeroState::dist_sample,
+            "sample particles for AeroState from an AeroDist")
     ;
 
     py::class_<GasData, std::shared_ptr<GasData>>(m, "GasData",

tests/test_aero_dist.py

@@ -13,12 +13,16 @@
 import PyPartMC as ppmc
 
 from .test_aero_data import AERO_DATA_CTOR_ARG_MINIMAL
-from .test_aero_mode import AERO_MODE_CTOR_LOG_NORMAL
+from .test_aero_mode import AERO_MODE_CTOR_LOG_NORMAL, AERO_MODE_CTOR_LOG_NORMAL_FULL
 
 AERO_DIST_CTOR_ARG_MINIMAL = [
     AERO_MODE_CTOR_LOG_NORMAL,
 ]
 
+AERO_DIST_CTOR_ARG_FULL = [
+    AERO_MODE_CTOR_LOG_NORMAL_FULL,
+]
+
 
 @pytest.fixture
 def sut_minimal():

tests/test_aero_mode.py

@@ -23,6 +23,17 @@
     }
 }
 
+AERO_MODE_CTOR_LOG_NORMAL_FULL = {
+    "test_mode": {
+        "mass_frac": [{"SO4": [1]}],
+        "diam_type": "geometric",
+        "mode_type": "log_normal",
+        "num_conc": 100 / si.m**3,
+        "geom_mean_diam": 2 * si.um,
+        "log10_geom_std_dev": np.log10(1.6),
+    }
+}
+
 
 class TestAeroMode:
     @staticmethod

tests/test_aero_state.py

@@ -13,15 +13,29 @@
 import PyPartMC as ppmc
 
 from .test_aero_data import AERO_DATA_CTOR_ARG_FULL, AERO_DATA_CTOR_ARG_MINIMAL
+from .test_aero_dist import AERO_DIST_CTOR_ARG_FULL, AERO_DIST_CTOR_ARG_MINIMAL
 from .test_env_state import ENV_STATE_CTOR_ARG_MINIMAL
 
 AERO_STATE_CTOR_ARG_MINIMAL = 44
 
 
 @pytest.fixture
 def sut_minimal():
+    aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+    aero_dist = ppmc.AeroDist(aero_data, AERO_DIST_CTOR_ARG_MINIMAL)
+    sut = ppmc.AeroState(AERO_STATE_CTOR_ARG_MINIMAL, aero_data)
+    _ = sut.dist_sample(aero_dist, 1.0, 0.0, True, True)
+    aero_data = None
+    gc.collect()
+    return sut
+
+
+@pytest.fixture
+def sut_full():
     aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_FULL)
+    aero_dist = ppmc.AeroDist(aero_data, AERO_DIST_CTOR_ARG_FULL)
     sut = ppmc.AeroState(AERO_STATE_CTOR_ARG_MINIMAL, aero_data)
+    _ = sut.dist_sample(aero_dist, 1.0, 0.0, True, True)
     aero_data = None
     gc.collect()
     return sut
@@ -40,18 +54,20 @@ def test_ctor():
         assert sut is not None
 
     @staticmethod
-    @pytest.mark.xfail(strict=True)  # TODO #116
-    @pytest.mark.parametrize("n_part", (1, 44, 666))
+    @pytest.mark.parametrize("n_part", (44, 666))
     def test_len(n_part):
         # arrange
         aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        aero_dist = ppmc.AeroDist(aero_data, AERO_DIST_CTOR_ARG_MINIMAL)
         sut = ppmc.AeroState(n_part, aero_data)
+        _ = sut.dist_sample(aero_dist, 1.0, 0.0, True, True)
 
         # act
         size = len(sut)
 
         # assert
-        assert size == n_part
+        assert int(size) > n_part * 0.5
+        assert int(size) < n_part * 2
 
     @staticmethod
     def test_copy(sut_minimal):  # pylint: disable=redefined-outer-name
@@ -128,18 +144,18 @@ def test_diameters(sut_minimal):  # pylint: disable=redefined-outer-name
         assert len(diameters) == len(sut_minimal)
 
     @staticmethod
-    def test_crit_rel_humids(sut_minimal):  # pylint: disable=redefined-outer-name
+    def test_crit_rel_humids(sut_full):  # pylint: disable=redefined-outer-name
         # arrange
         args = {"rel_humidity": 0.8, **ENV_STATE_CTOR_ARG_MINIMAL}
         env_state = ppmc.EnvState(args)
         env_state.set_temperature(300)
 
         # act
-        crit_rel_humids = sut_minimal.crit_rel_humids(env_state)
+        crit_rel_humids = sut_full.crit_rel_humids(env_state)
 
         # assert
         assert isinstance(crit_rel_humids, list)
-        assert len(crit_rel_humids) == len(sut_minimal)
+        assert len(crit_rel_humids) == len(sut_full)
         assert (np.asarray(crit_rel_humids) > 1).all()
         assert (np.asarray(crit_rel_humids) < 1.2).all()
 
@@ -155,7 +171,7 @@ def test_mixing_state(sut_minimal):  # pylint: disable=redefined-outer-name
     @staticmethod
     def test_bin_average_comp(sut_minimal):  # pylint: disable=redefined-outer-name
         # arrange
-        bin_grid = ppmc.BinGrid(123, "log", 1, 100)
+        bin_grid = ppmc.BinGrid(123, "log", 1e-9, 1e-4)
 
         # act
         sut_minimal.bin_average_comp(bin_grid)
@@ -214,3 +230,14 @@ def test_get_particle_out_of_range(
 
         # assert
         assert False
+
+    @staticmethod
+    def test_dist_sample():
+        n_part = 44
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        aero_dist = ppmc.AeroDist(aero_data, AERO_DIST_CTOR_ARG_MINIMAL)
+        sut = ppmc.AeroState(n_part, aero_data)
+        n_added = sut.dist_sample(aero_dist, 1.0, 0.0, True, True)
+
+        assert n_added > n_part * 0.5
+        assert n_added < n_part * 2