From fa6007f96259e49c2ec57f5dd7fb7a18b25f5c75 Mon Sep 17 00:00:00 2001 From: Zach D'Aquino <107204919+zdaq12@users.noreply.github.com> Date: Wed, 17 May 2023 15:20:43 -0500 Subject: [PATCH] Aero particle additions (#212) Co-authored-by: Sylwester Arabas --- src/aero_particle.F90 | 196 +++++++++++++++++++++++++++++++++- src/aero_particle.hpp | 106 +++++++++++++++++- src/pypartmc.cpp | 18 ++++ tests/test_aero_particle.py | 207 +++++++++++++++++++++++++++++++++++- 4 files changed, 523 insertions(+), 4 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 5d67c19e..41abf5b6 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -36,7 +36,7 @@ subroutine f_aero_particle_init(ptr_c, aero_data_ptr_c, arr_data, arr_size) bind real(c_double), dimension(arr_size), intent(in) :: arr_data call c_f_pointer(ptr_c, ptr_f) call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) - allocate(ptr_f%vol(arr_size)) + call aero_particle_zero(ptr_f, aero_data_ptr_f) call aero_particle_set_vols(ptr_f, arr_data) end subroutine @@ -193,4 +193,198 @@ subroutine f_aero_particle_solute_kappa(aero_particle_ptr_c, aero_data_ptr_c, ka kappa = aero_particle_solute_kappa(aero_particle_ptr_f, aero_data_ptr_f) end subroutine + subroutine f_aero_particle_moles(aero_particle_ptr_c, aero_data_ptr_c, moles) bind(C) + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(aero_data_t), pointer :: aero_data_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c, aero_data_ptr_c + real(c_double), intent(out) :: moles + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) + + moles = aero_particle_moles(aero_particle_ptr_f, aero_data_ptr_f) + end subroutine + + subroutine f_aero_particle_mobility_diameter( & + aero_particle_ptr_c, & + aero_data_ptr_c, & + env_state_ptr_c, & + mobility_diameter & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(aero_data_t), pointer :: aero_data_ptr_f => null() + type(env_state_t), pointer :: env_state_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c, aero_data_ptr_c, & + env_state_ptr_c + real(c_double), intent(out) :: mobility_diameter + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) + call c_f_pointer(env_state_ptr_c, env_state_ptr_f) + + mobility_diameter = aero_particle_mobility_diameter( & + aero_particle_ptr_f, & + aero_data_ptr_f, & + env_state_ptr_f & + ) + end subroutine + + subroutine f_aero_particle_density( & + aero_particle_ptr_c, & + aero_data_ptr_c, & + density & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(aero_data_t), pointer :: aero_data_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c, aero_data_ptr_c + real(c_double), intent(out) :: density + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) + + density = aero_particle_density( & + aero_particle_ptr_f, & + aero_data_ptr_f & + ) + end subroutine + + subroutine f_aero_particle_approx_crit_rel_humid( & + aero_particle_ptr_c, & + aero_data_ptr_c, & + env_state_ptr_c, & + approx_crit_rel_humid & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(aero_data_t), pointer :: aero_data_ptr_f => null() + type(env_state_t), pointer :: env_state_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c, aero_data_ptr_c, & + env_state_ptr_c + real(c_double), intent(out) :: approx_crit_rel_humid + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) + call c_f_pointer(env_state_ptr_c, env_state_ptr_f) + + approx_crit_rel_humid = aero_particle_approx_crit_rel_humid( & + aero_particle_ptr_f, & + aero_data_ptr_f, & + env_state_ptr_f & + ) + end subroutine + + subroutine f_aero_particle_crit_rel_humid( & + aero_particle_ptr_c, & + aero_data_ptr_c, & + env_state_ptr_c, & + crit_rel_humid & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(aero_data_t), pointer :: aero_data_ptr_f => null() + type(env_state_t), pointer :: env_state_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c, aero_data_ptr_c, & + env_state_ptr_c + real(c_double), intent(out) :: crit_rel_humid + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) + call c_f_pointer(env_state_ptr_c, env_state_ptr_f) + + crit_rel_humid = aero_particle_crit_rel_humid( & + aero_particle_ptr_f, & + aero_data_ptr_f, & + env_state_ptr_f & + ) + end subroutine + + subroutine f_aero_particle_crit_diameter( & + aero_particle_ptr_c, & + aero_data_ptr_c, & + env_state_ptr_c, & + crit_diameter & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(aero_data_t), pointer :: aero_data_ptr_f => null() + type(env_state_t), pointer :: env_state_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c, aero_data_ptr_c, & + env_state_ptr_c + real(c_double), intent(out) :: crit_diameter + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) + call c_f_pointer(env_state_ptr_c, env_state_ptr_f) + + crit_diameter = aero_particle_crit_diameter( & + aero_particle_ptr_f, & + aero_data_ptr_f, & + env_state_ptr_f & + ) + end subroutine + + subroutine f_aero_particle_coagulate( & + aero_particle_1_ptr_c, & + aero_particle_2_ptr_c, & + aero_particle_new_ptr_c & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_1_ptr_f => null() + type(aero_particle_t), pointer :: aero_particle_2_ptr_f => null() + type(aero_particle_t), pointer :: aero_particle_new_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_1_ptr_c, aero_particle_2_ptr_c + type(c_ptr), intent(inout) :: aero_particle_new_ptr_c + + call c_f_pointer(aero_particle_1_ptr_c, aero_particle_1_ptr_f) + call c_f_pointer(aero_particle_2_ptr_c, aero_particle_2_ptr_f) + call c_f_pointer(aero_particle_new_ptr_c, aero_particle_new_ptr_f) + + call aero_particle_coagulate( & + aero_particle_1_ptr_f, & + aero_particle_2_ptr_f, & + aero_particle_new_ptr_f & + ) + + aero_particle_new_ptr_c = c_loc(aero_particle_new_ptr_f) + end subroutine + + subroutine f_aero_particle_zero( & + aero_particle_ptr_c, & + aero_data_ptr_c & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(aero_data_t), pointer :: aero_data_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c, aero_data_ptr_c + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f) + + call aero_particle_zero( & + aero_particle_ptr_f, & + aero_data_ptr_f & + ) + end subroutine + + subroutine f_aero_particle_set_vols( & + aero_particle_ptr_c, & + vol_size, & + volumes & + ) bind(C) + + type(aero_particle_t), pointer :: aero_particle_ptr_f => null() + type(c_ptr), intent(in) :: aero_particle_ptr_c + integer(c_int), intent(in) :: vol_size + real(c_double), dimension(vol_size), intent(in) :: volumes + + call c_f_pointer(aero_particle_ptr_c, aero_particle_ptr_f) + + call aero_particle_set_vols( & + aero_particle_ptr_f, & + volumes & + ) + end subroutine + end module diff --git a/src/aero_particle.hpp b/src/aero_particle.hpp index 2924fed2..1250d957 100644 --- a/src/aero_particle.hpp +++ b/src/aero_particle.hpp @@ -8,6 +8,7 @@ #include "pmc_resource.hpp" #include "aero_data.hpp" +#include "env_state.hpp" #include "pybind11/stl.h" extern "C" void f_aero_particle_ctor(void *ptr) noexcept; @@ -25,7 +26,15 @@ extern "C" void f_aero_particle_mass(const void *aero_particle_ptr, const void * extern "C" void f_aero_particle_species_mass(const void *aero_particle_ptr, const int *i_spec, const void *aero_data_ptr, double *mass) noexcept; extern "C" void f_aero_particle_species_masses(const void *aero_particle_ptr, const void *aero_data_ptr, const int *size_masses, void *masses) noexcept; extern "C" void f_aero_particle_solute_kappa(const void *aero_particle_ptr, const void *aero_data_ptr, void *kappa) noexcept; - +extern "C" void f_aero_particle_moles(const void *aero_particle_ptr, const void *aero_data_ptr, void *moles) noexcept; +extern "C" void f_aero_particle_mobility_diameter(const void *aero_particle_ptr, const void *aero_data_ptr, const void *env_state_ptr, void *mobility_diameter) noexcept; +extern "C" void f_aero_particle_density(const void *aero_particle_ptr, const void *aero_data_ptr, void *density) noexcept; +extern "C" void f_aero_particle_approx_crit_rel_humid(const void *aero_particle_ptr, const void *aero_data_ptr, const void *env_state_ptr, void *approx_crit_rel_humid) noexcept; +extern "C" void f_aero_particle_crit_rel_humid(const void *aero_particle_ptr, const void *aero_data_ptr, const void *env_state_ptr, void *crit_rel_humid) noexcept; +extern "C" void f_aero_particle_crit_diameter(const void *aero_particle_ptr, const void *aero_data_ptr, const void *env_state, void *crit_diameter) noexcept; +extern "C" void f_aero_particle_coagulate(const void *aero_particle_1_ptr, const void *aero_particle_2_ptr, void *new_particle_ptr) noexcept; +extern "C" void f_aero_particle_zero(void *aero_particle_ptr, const void *aero_data_ptr) noexcept; +extern "C" void f_aero_particle_set_vols(void *aero_particle_ptr, const int *vol_size, const void *volumes) noexcept; namespace py = pybind11; struct AeroParticle { @@ -169,5 +178,98 @@ struct AeroParticle { return kappa; } -}; + static auto moles(const AeroParticle &self) { + double moles; + f_aero_particle_moles( + self.ptr.f_arg(), + self.aero_data.get(), + &moles + ); + return moles; + } + + static auto mobility_diameter(const AeroParticle &self, const EnvState &env_state) { + double mobility_diameter; + f_aero_particle_mobility_diameter( + self.ptr.f_arg(), + self.aero_data.get(), + env_state.ptr.f_arg(), + &mobility_diameter + ); + return mobility_diameter; + } + + static auto density(const AeroParticle &self) { + double density; + f_aero_particle_density( + self.ptr.f_arg(), + self.aero_data.get(), + &density + ); + return density; + } + + static auto approx_crit_rel_humid(const AeroParticle &self, const EnvState &env_state) { + double approx_crit_rel_humid; + f_aero_particle_approx_crit_rel_humid( + self.ptr.f_arg(), + self.aero_data.get(), + env_state.ptr.f_arg(), + &approx_crit_rel_humid + ); + return approx_crit_rel_humid; + } + + static auto crit_rel_humid(const AeroParticle &self, const EnvState &env_state) { + double crit_rel_humid; + f_aero_particle_crit_rel_humid( + self.ptr.f_arg(), + self.aero_data.get(), + env_state.ptr.f_arg(), + &crit_rel_humid + ); + return crit_rel_humid; + } + + static auto crit_diameter(const AeroParticle &self, const EnvState &env_state) { + double crit_diameter; + f_aero_particle_crit_diameter( + self.ptr.f_arg(), + self.aero_data.get(), + env_state.ptr.f_arg(), + &crit_diameter + ); + return crit_diameter; + } + + static auto coagulate(const AeroParticle &self, const AeroParticle &two) { + int len = AeroData::__len__(*self.aero_data); + std::valarray data(len); + AeroParticle* new_ptr = new AeroParticle(self.aero_data, data); + f_aero_particle_coagulate( + self.ptr.f_arg(), + two.ptr.f_arg(), + new_ptr + ); + return new_ptr; + } + + static void zero(AeroParticle &self) { + f_aero_particle_zero( + self.ptr.f_arg_non_const(), + self.aero_data.get() + ); + } + + static void set_vols(AeroParticle &self, const std::valarray&volumes) { + int len = AeroData::__len__(*self.aero_data.get()); + if (volumes.size() != size_t(len)) + throw std::runtime_error("AeroData size mistmatch"); + f_aero_particle_set_vols( + self.ptr.f_arg_non_const(), + &len, + begin(volumes) + ); + } +}; \ No newline at end of file diff --git a/src/pypartmc.cpp b/src/pypartmc.cpp index 1f0e53e4..ead4dc82 100644 --- a/src/pypartmc.cpp +++ b/src/pypartmc.cpp @@ -124,6 +124,24 @@ PYBIND11_MODULE(_PyPartMC, m) { "Mass of all species in the particle (kg).") .def_property_readonly("solute_kappa", AeroParticle::solute_kappa, "Returns the average of the solute kappas (1).") + .def_property_readonly("moles", AeroParticle::moles, + "Total moles in the particle (1).") + .def("mobility_diameter", AeroParticle::mobility_diameter, + "Mobility diameter of the particle (m).") + .def_property_readonly("density", AeroParticle::density, + "Average density of the particle (kg/m^3)") + .def("approx_crit_rel_humid", AeroParticle::approx_crit_rel_humid, + "Returns the approximate critical relative humidity (1).") + .def("crit_rel_humid", AeroParticle::crit_rel_humid, + "Returns the critical relative humidity (1).") + .def("crit_diameter", AeroParticle::crit_diameter, + "Returns the critical diameter (m).") + .def("coagulate", AeroParticle::coagulate, + "Coagulate two particles together to make a new one. The new particle will not have its ID set.") + .def("zero", AeroParticle::zero, + "Resets an aero_particle to be zero.") + .def("set_vols", AeroParticle::set_vols, + "Sets the aerosol particle volumes.") ; py::class_(m, "AeroState", diff --git a/tests/test_aero_particle.py b/tests/test_aero_particle.py index 21e9e1ac..17e83913 100644 --- a/tests/test_aero_particle.py +++ b/tests/test_aero_particle.py @@ -13,9 +13,10 @@ from PyPartMC import si from .test_aero_data import AERO_DATA_CTOR_ARG_MINIMAL +from .test_env_state import ENV_STATE_CTOR_ARG_MINIMAL -class TestAeroParticle: +class TestAeroParticle: # pylint: disable=too-many-public-methods @staticmethod @pytest.mark.parametrize( "volumes", @@ -291,3 +292,207 @@ def test_solute_kappa(): # assert np.testing.assert_almost_equal(kappa, 1.479240661) + + @staticmethod + def test_moles(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + aero_data = None + volumes = None + gc.collect() + + # act + moles = sut.moles + check = 1000 / 18e-3 + 4400 / 35.5e-3 + 6600 / 23e-3 + + # assert + np.testing.assert_almost_equal(moles, check) + + @staticmethod + def test_mobility_diameter(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + env_state = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL) + aero_data = None + volumes = None + + # act + mobility_diameter = sut.mobility_diameter(env_state) + env_state = None + gc.collect() + + # assert + assert mobility_diameter is not None + + @staticmethod + def test_density(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + aero_data = None + volumes = None + gc.collect() + + # act + density = sut.density + + # assert + assert density == (1000 * 1 + 2200 * 2 + 2200 * 3) / 6 + + @staticmethod + def test_approx_crit_rel_humid(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + env_state = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL) + env_state.set_temperature(288) + aero_data = None + volumes = None + + # act + approx_crit_rel_humid = sut.approx_crit_rel_humid(env_state) + env_state = None + gc.collect() + + # assert + assert approx_crit_rel_humid is not None + + @staticmethod + def test_crit_rel_humid(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + env_state = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL) + env_state.set_temperature(288) + aero_data = None + volumes = None + + # act + crit_rel_humid = sut.crit_rel_humid(env_state) + env_state = None + gc.collect() + + # assert + assert crit_rel_humid is not None + + @staticmethod + def test_crit_diameter(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + env_state = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL) + env_state.set_temperature(288) + aero_data = None + volumes = None + + # act + crit_diameter = sut.crit_diameter(env_state) + env_state = None + gc.collect() + + # assert + assert crit_diameter is not None + + @staticmethod + def test_coagulate(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes_1 = [1, 2, 3] + volumes_2 = [3, 2, 1] + sut = ppmc.AeroParticle(aero_data, volumes_1) + aero_particle_2 = ppmc.AeroParticle(aero_data, volumes_2) + aero_data = None + volumes_1 = None + volumes_2 = None + gc.collect() + + # act + coagulated = sut.coagulate(aero_particle_2) + + # assert + assert coagulated.volumes == [4, 4, 4] + + @staticmethod + def test_zero(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + aero_data = None + volumes = None + gc.collect() + + # act + sut.zero() + + # assert + assert sut.volumes == [0, 0, 0] + + @staticmethod + def test_set_vols(): + # arrange + aero_data_arg = ( + {"H2O": [1000 * si.kg / si.m**3, 0, 18e-3 * si.kg / si.mol, 0]}, + {"Cl": [2200 * si.kg / si.m**3, 1, 35.5e-3 * si.kg / si.mol, 0]}, + {"Na": [2200 * si.kg / si.m**3, 1, 23e-3 * si.kg / si.mol, 0]}, + ) + aero_data = ppmc.AeroData(aero_data_arg) + volumes = [1, 2, 3] + sut = ppmc.AeroParticle(aero_data, volumes) + aero_data = None + volumes = None + gc.collect() + + # act + sut.set_vols([3, 2, 1]) + + # assert + assert sut.volumes == [3, 2, 1]