Aero particle additions (#212)

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

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

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<double> 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<double>&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)
+ );
+ }
 
+};

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_<AeroState>(m, "AeroState",