From d6dbe31604174f56f5772f864af229b307e86d4b Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 22 May 2024 23:07:45 +0200 Subject: [PATCH] support for AeroMode instantiation with sampled mode type (incl. fixes in spec file logic) + sampled and mono mode types coverage in unit tests (in test_aero_[mode,state,dist].py) (#357) Co-authored-by: Jeffrey Curtis --- src/aero_mode.F90 | 36 ++++++++++ src/aero_mode.hpp | 59 +++++++++++++++- src/json_resource.hpp | 12 ++-- src/pypartmc.cpp | 4 ++ src/spec_file_pypartmc.F90 | 2 + src/spec_file_pypartmc.cpp | 10 +-- tests/test_aero_dist.py | 19 ++++++ tests/test_aero_mode.py | 136 +++++++++++++++++++++++++++++++++++++ tests/test_aero_state.py | 48 +++++++++++++ 9 files changed, 316 insertions(+), 10 deletions(-) diff --git a/src/aero_mode.F90 b/src/aero_mode.F90 index 1b514205..91952c10 100644 --- a/src/aero_mode.F90 +++ b/src/aero_mode.F90 @@ -244,4 +244,40 @@ subroutine f_aero_mode_get_name(ptr_c, name_data, name_size) bind(C) name_data = c_loc(aero_mode%name) name_size = len_trim(aero_mode%name) end subroutine + + subroutine f_aero_mode_get_sample_num_conc(ptr_c, arr_data, data_size) bind(C) + type(c_ptr), intent(in) :: ptr_c + type(aero_mode_t), pointer :: aero_mode => null() + integer(c_int), intent(in) :: data_size + real(c_double), dimension(data_size), intent(inout) :: arr_data + + call c_f_pointer(ptr_c, aero_mode) + + arr_data = aero_mode%sample_num_conc + + end subroutine + + subroutine f_aero_mode_get_sample_radius(ptr_c, arr_data, data_size) bind(C) + type(c_ptr), intent(in) :: ptr_c + type(aero_mode_t), pointer :: aero_mode => null() + integer(c_int), intent(in) :: data_size + real(c_double), dimension(data_size), intent(inout) :: arr_data + + call c_f_pointer(ptr_c, aero_mode) + + arr_data = aero_mode%sample_radius + + end subroutine + + subroutine f_aero_mode_get_sample_bins(ptr_c, n_bins) bind(c) + type(c_ptr), intent(in) :: ptr_c + type(aero_mode_t), pointer :: aero_mode => null() + integer(c_int), intent(out) :: n_bins + + call c_f_pointer(ptr_c, aero_mode) + + n_bins = size(aero_mode%sample_num_conc) + + end subroutine + end module diff --git a/src/aero_mode.hpp b/src/aero_mode.hpp index d2e4685a..e7576d87 100644 --- a/src/aero_mode.hpp +++ b/src/aero_mode.hpp @@ -118,6 +118,20 @@ extern "C" void f_aero_mode_from_json( void *aero_data_ptr ) noexcept; +extern "C" void f_aero_mode_get_sample_num_conc( + const void *ptr, + void *sample_num_conc_data, + const int *sample_num_conc_data_size +) noexcept; + +extern "C" void f_aero_mode_get_sample_radius( + const void *ptr, + void *sample_radius_data, + const int *sample_radius_data_size +) noexcept; + +extern "C" void f_aero_mode_get_sample_bins(const void *ptr, int *n_bins) noexcept; + struct AeroMode { PMCResource ptr; @@ -145,6 +159,24 @@ struct AeroMode { throw std::runtime_error("mass_frac value must be a list of single-element dicts"); if (!InputJSONResource::unique_keys(mass_frac)) throw std::runtime_error("mass_frac keys must be unique"); + if (mode["mode_type"] == "sampled") { + if (mode.find("size_dist") == mode.end()) + throw std::runtime_error("size_dist key must be set for mode_type=sampled"); + auto sd = mode["size_dist"]; + if ( + sd.size() != 2 || + !sd[0].is_object() || + sd[0].size() != 1 || + sd[1].size() != 1 || + sd[0].find("diam") == sd[0].end() || + sd[1].find("num_conc") == sd[1].end() + ) + throw std::runtime_error("size_dist value must be an iterable of two single-element dicts (first with 'diam', second with 'num_conc' as keys)"); + auto diam = *sd[0].find("diam"); + auto num_conc = *sd[1].find("num_conc"); + if (diam.size() != num_conc.size() + 1) + throw std::runtime_error("size_dist['num_conc'] must have len(size_dist['diam'])-1 elements"); + } } static auto get_num_conc(const AeroMode &self){ @@ -291,7 +323,7 @@ struct AeroMode { int type; f_aero_mode_get_type(self.ptr.f_arg(), &type); - if (type < 0 || (unsigned int)type >= AeroMode::types().size()) + if (type <= 0 || (unsigned int)type > AeroMode::types().size()) throw std::logic_error("Unknown mode type."); return AeroMode::types()[type - 1]; @@ -312,4 +344,29 @@ struct AeroMode { name[i] = f_ptr[i]; return name; } + + static auto get_sample_radius(const AeroMode &self) { + int len; + f_aero_mode_get_sample_bins(self.ptr.f_arg(), &len); + len++; + std::valarray sample_radius(len); + f_aero_mode_get_sample_radius( + self.ptr.f_arg(), + begin(sample_radius), + &len + ); + return sample_radius; + } + + static auto get_sample_num_conc(const AeroMode &self) { + int len; + f_aero_mode_get_sample_bins(self.ptr.f_arg(), &len); + std::valarray sample_num_conc(len); + f_aero_mode_get_sample_num_conc( + self.ptr.f_arg(), + begin(sample_num_conc), + &len + ); + return sample_num_conc; + } }; diff --git a/src/json_resource.hpp b/src/json_resource.hpp index e3736b7e..91e46219 100644 --- a/src/json_resource.hpp +++ b/src/json_resource.hpp @@ -26,7 +26,7 @@ struct JSONResource { } protected: - size_t index = 0; + size_t index = 0, named_array_read_count = 0; JSONResource() {} @@ -71,6 +71,10 @@ struct JSONResource { public: virtual ~JSONResource() {} + auto n_named_array_read_count() noexcept { + return this->named_array_read_count; + } + void zoom_in(const bpstd::string_view &sub) noexcept { auto it = this->json->is_array() ? this->json->at(this->json->size()-1).find(sub) @@ -84,6 +88,7 @@ struct JSONResource { else this->set_current_json_ptr(&(*it)); + this->named_array_read_count = 0; } void zoom_out() noexcept { @@ -101,13 +106,12 @@ struct JSONResource { return this->json->begin(); } - // TODO #112: to be removed after initialising GasData with a list, and not JSON? - auto first_field_name() const noexcept { + auto first_field_name() noexcept { // TODO #112: handle errors std::string name = ""; assert(this->json->size() > 0); assert(this->json->begin()->size() > 0); - for (auto &entry : this->json->at(0).items()) + for (auto &entry : this->json->at(this->named_array_read_count++).items()) { name = entry.key(); } diff --git a/src/pypartmc.cpp b/src/pypartmc.cpp index 900b4726..94b200b6 100644 --- a/src/pypartmc.cpp +++ b/src/pypartmc.cpp @@ -453,6 +453,10 @@ PYBIND11_MODULE(_PyPartMC, m) { .def_property("gsd", &AeroMode::get_gsd, &AeroMode::set_gsd, "Geometric standard deviation") .def("set_sample", &AeroMode::set_sampled) + .def_property_readonly("sample_num_conc", &AeroMode::get_sample_num_conc, + "Sample bin number concentrations (m^{-3})") + .def_property_readonly("sample_radius", &AeroMode::get_sample_radius, + "Sample bin radii (m).") .def_property("type", &AeroMode::get_type, &AeroMode::set_type, "Mode type (given by module constants)") .def_property("name", &AeroMode::get_name, &AeroMode::set_name, diff --git a/src/spec_file_pypartmc.F90 b/src/spec_file_pypartmc.F90 index 31116516..a92e7433 100644 --- a/src/spec_file_pypartmc.F90 +++ b/src/spec_file_pypartmc.F90 @@ -129,6 +129,8 @@ subroutine spec_file_read_real_named_array(file, max_lines, names, vals) n_rows = max_lines end if + if (allocated(names)) deallocate(names) + if (allocated(vals)) deallocate(vals) allocate(names(n_rows)) allocate(vals(n_rows, n_cols)) allocate(vals_row(max(1, n_cols))) diff --git a/src/spec_file_pypartmc.cpp b/src/spec_file_pypartmc.cpp index fdda521f..214d94f0 100644 --- a/src/spec_file_pypartmc.cpp +++ b/src/spec_file_pypartmc.cpp @@ -145,10 +145,10 @@ void spec_file_read_real_named_array_size( int *n_rows, int *n_cols ) noexcept { - auto first_field = json_resource_ptr()->first_field_name(); *n_rows = json_resource_ptr()->n_numeric_array_entries(); + + auto first_field = json_resource_ptr()->first_field_name(); *n_cols = json_resource_ptr()->n_elements(first_field); - // TODO #112: check each line has the same number of elements as time } extern "C" @@ -174,17 +174,17 @@ void spec_file_read_real_named_array_data( ++i, ++it ) { assert(it->is_object()); - if (i == row-1) { + if (i == (row - 1) + (json_resource_ptr()->n_named_array_read_count() - 1)) { assert(it->size() == 1); for (auto &entry : it->items()) { - // TODO #112: use input name_size as limit param + assert(*name_size > (long)entry.key().size()); for (auto c=0u; c < entry.key().size(); ++c) name_data[c] = entry.key()[c]; *name_size = entry.key().size(); for (auto idx=0u; idx < entry.value().size(); ++idx) { vals[idx] = entry.value().at(idx).get(); } - break; // TODO #112 + break; } } } diff --git a/tests/test_aero_dist.py b/tests/test_aero_dist.py index adb64cac..94e2d6eb 100644 --- a/tests/test_aero_dist.py +++ b/tests/test_aero_dist.py @@ -217,3 +217,22 @@ def test_ctor_error_on_repeated_massfrac_keys(): # assert assert str(exc_info.value) == "mass_frac keys must be unique" + + @staticmethod + def test_ctor_sampled_mode(): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + ctor_arg = copy.deepcopy(AERO_DIST_CTOR_ARG_MINIMAL) + ctor_arg[0]["test_mode"]["mode_type"] = "sampled" + ctor_arg[0]["test_mode"]["size_dist"] = [ + {"diam": [1, 2, 3, 4]}, + {"num_conc": [1, 2, 3]}, + ] + + # act + sut = ppmc.AeroDist(aero_data, ctor_arg) + + # assert + assert sut.mode(0).num_conc == sum( + ctor_arg[0]["test_mode"]["size_dist"][1]["num_conc"] + ) diff --git a/tests/test_aero_mode.py b/tests/test_aero_mode.py index fb148698..1e0b1c18 100644 --- a/tests/test_aero_mode.py +++ b/tests/test_aero_mode.py @@ -48,6 +48,18 @@ } } +AERO_MODE_CTOR_SAMPLED = { + "test_mode": { + "mass_frac": [{"H2O": [1]}], + "diam_type": "geometric", + "mode_type": "sampled", + "size_dist": [ + {"diam": [1, 2, 3, 4]}, + {"num_conc": [100, 200, 300]}, + ], + } +} + class TestAeroMode: @staticmethod @@ -289,3 +301,127 @@ def test_segfault_case(): # TODO #319 ) print(fishy_ctor_arg) ppmc.AeroMode(aero_data, fishy_ctor_arg) + + @staticmethod + @pytest.mark.skipif(platform.machine() == "arm64", reason="TODO #348") + def test_sampled_without_size_dist(): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + fishy_ctor_arg = copy.deepcopy(AERO_MODE_CTOR_LOG_NORMAL) + fishy_ctor_arg["test_mode"]["mode_type"] = "sampled" + + # act + with pytest.raises(Exception) as exc_info: + ppmc.AeroMode(aero_data, fishy_ctor_arg) + + # assert + assert str(exc_info.value) == "size_dist key must be set for mode_type=sampled" + + @staticmethod + @pytest.mark.parametrize( + "fishy", + ( + None, + [], + [{}, {}, {}], + [{}, []], + [{"diam": None}, {}], + [{"num_conc": None}, {}], + [{"diam": None, "": None}, {}], + [{"num_conc": None, "": None}, {}], + ), + ) + @pytest.mark.skipif(platform.machine() == "arm64", reason="TODO #348") + def test_sampled_with_fishy_size_dist(fishy): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + fishy_ctor_arg = copy.deepcopy(AERO_MODE_CTOR_LOG_NORMAL) + fishy_ctor_arg["test_mode"]["mode_type"] = "sampled" + fishy_ctor_arg["test_mode"]["size_dist"] = fishy + + # act + with pytest.raises(Exception) as exc_info: + ppmc.AeroMode(aero_data, fishy_ctor_arg) + + # assert + assert ( + str(exc_info.value) + == "size_dist value must be an iterable of two single-element dicts" + + " (first with 'diam', second with 'num_conc' as keys)" + ) + + @staticmethod + @pytest.mark.skipif(platform.machine() == "arm64", reason="TODO #348") + def test_sampled_with_diam_of_different_len_than_num_conc(): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + fishy_ctor_arg = copy.deepcopy(AERO_MODE_CTOR_LOG_NORMAL) + fishy_ctor_arg["test_mode"]["mode_type"] = "sampled" + fishy_ctor_arg["test_mode"]["size_dist"] = [ + {"diam": [1, 2, 3]}, + {"num_conc": [1, 2, 3]}, + ] + + # act + with pytest.raises(Exception) as exc_info: + ppmc.AeroMode(aero_data, fishy_ctor_arg) + + # assert + assert ( + str(exc_info.value) + == "size_dist['num_conc'] must have len(size_dist['diam'])-1 elements" + ) + + @staticmethod + def test_sampled(): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + + # act + sut = ppmc.AeroMode(aero_data, AERO_MODE_CTOR_SAMPLED) + + # assert + assert sut.type == "sampled" + assert sut.num_conc == np.sum( + AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][1]["num_conc"] + ) + assert ( + sut.sample_num_conc + == AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][1]["num_conc"] + ) + assert ( + np.array(sut.sample_radius) * 2 + == AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][0]["diam"] + ).all() + + @staticmethod + def test_set_sample(): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + + diams = [1, 2, 3, 4] + num_concs = [100, 200, 300] + sut = ppmc.AeroMode( + aero_data, + { + "test_mode": { + "mass_frac": [{"H2O": [1]}], + "diam_type": "geometric", + "mode_type": "sampled", + "size_dist": [ + {"diam": diams}, + {"num_conc": num_concs}, + ], + } + }, + ) + num_conc_orig = sut.num_conc + # act + diams = [0.5 * x for x in diams] + num_concs = [2 * x for x in num_concs] + sut.set_sample(diams, num_concs) + # assert + assert sut.num_conc == np.sum(num_concs) + assert sut.sample_num_conc == num_concs + assert (np.array(sut.sample_radius) * 2 == diams).all() + assert sut.num_conc == num_conc_orig * 2 diff --git a/tests/test_aero_state.py b/tests/test_aero_state.py index 85217bd6..74813d36 100644 --- a/tests/test_aero_state.py +++ b/tests/test_aero_state.py @@ -18,6 +18,7 @@ AERO_DIST_CTOR_ARG_FULL, AERO_DIST_CTOR_ARG_MINIMAL, ) +from .test_aero_mode import AERO_MODE_CTOR_SAMPLED from .test_env_state import ENV_STATE_CTOR_ARG_MINIMAL AERO_STATE_CTOR_ARG_MINIMAL = 44, "nummass_source" @@ -508,3 +509,50 @@ def test_dist_sample(args, weighting): # assert assert n_added > n_part * 0.5 assert n_added < n_part * 2 + + @staticmethod + def test_dist_sample_sampled(): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + aero_dist = ppmc.AeroDist(aero_data, [AERO_MODE_CTOR_SAMPLED]) + sut = ppmc.AeroState(aero_data, *AERO_STATE_CTOR_ARG_MINIMAL) + + # act + _ = sut.dist_sample(aero_dist, 1.0, 0.0, True, True) + + # assert + assert ( + np.array(sut.diameters()) + >= AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][0]["diam"][0] + ).all() + assert ( + np.array(sut.diameters()) + <= AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][0]["diam"][-1] + ).all() + + @staticmethod + def test_dist_sample_mono(): + # arrange + aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL) + diam = 2e-6 + aero_dist = ppmc.AeroDist( + aero_data, + [ + { + "test_mode": { + "mass_frac": [{"H2O": [1]}], + "diam_type": "geometric", + "mode_type": "mono", + "num_conc": 1e12, + "diam": diam, + } + } + ], + ) + sut = ppmc.AeroState(aero_data, *AERO_STATE_CTOR_ARG_MINIMAL) + + # act + _ = sut.dist_sample(aero_dist, 1.0, 0.0, True, True) + + # assert + assert np.isclose(np.array(sut.diameters()), diam).all()