Merge pull request #143 from zdaq12/dry_dep

Add loss function support

Author: slayoo

src/env_state.F90

@@ -69,4 +69,48 @@ subroutine f_env_state_get_rel_humid(ptr_c, rel_humid) bind(C)
 
     end subroutine
 
+    subroutine f_env_state_set_height(ptr_c, height) bind(C)
+        type(env_state_t), pointer :: ptr_f => null()
+        type(c_ptr), intent(in) :: ptr_c
+        real(c_double), intent(in) :: height
+
+        call c_f_pointer(ptr_c, ptr_f)
+
+        ptr_f%height = height
+
+    end subroutine
+
+    subroutine f_env_state_get_height(ptr_c, height) bind(C)
+        type(env_state_t), pointer :: ptr_f => null()
+        type(c_ptr), intent(in) :: ptr_c
+        real(c_double), intent(out) :: height
+        
+        call c_f_pointer(ptr_c, ptr_f)
+
+        height = ptr_f%height
+
+    end subroutine
+
+    subroutine f_env_state_set_pressure(ptr_c, pressure) bind(C)
+        type(env_state_t), pointer :: ptr_f => null()
+        type(c_ptr), intent(in) :: ptr_c
+        real(c_double), intent(in) :: pressure
+
+        call c_f_pointer(ptr_c, ptr_f)
+
+        ptr_f%pressure = pressure
+
+    end subroutine
+
+    subroutine f_env_state_get_pressure(ptr_c, pressure) bind(C)
+        type(env_state_t), pointer :: ptr_f => null()
+        type(c_ptr), intent(in) :: ptr_c
+        real(c_double), intent(out) :: pressure
+
+        call c_f_pointer(ptr_c, ptr_f)
+
+        pressure = ptr_f%pressure
+
+    end subroutine
+
 end module

src/env_state.hpp

@@ -15,6 +15,10 @@ extern "C" void f_env_state_from_json(const void *ptr) noexcept;
 extern "C" void f_env_state_set_temperature(const void *ptr, const double *temperature) noexcept;
 extern "C" void f_env_state_get_temperature(const void *ptr, double *temperature) noexcept;
 extern "C" void f_env_state_get_rel_humid(const void *ptr, double *rel_humid) noexcept;
+extern "C" void f_env_state_set_height(const void *ptr, const double *height) noexcept;
+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;
 
 struct EnvState {
     PMCResource ptr;
@@ -44,4 +48,26 @@ struct EnvState {
         f_env_state_get_rel_humid(&self.ptr, &rel_humid);
         return rel_humid;
     }
+    
+    static void set_height(const EnvState &self, const double height) {
+        f_env_state_set_height(&self.ptr, &height);
+    }
+
+    static double get_height(const EnvState &self) {
+        double height;
+
+        f_env_state_get_height(&self.ptr, &height);
+        return height;
+    }
+
+    static void set_pressure(const EnvState &self, const double pressure) {
+        f_env_state_set_pressure(&self.ptr, &pressure);
+    }
+
+    static double get_pressure(const EnvState &self) {
+        double pressure;
+
+        f_env_state_get_pressure(&self.ptr, &pressure);
+        return pressure;
+    }
 };

src/pypartmc.cpp

@@ -168,6 +168,8 @@ 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("height", &EnvState::get_height, &EnvState::set_height)
+        .def_property("pressure", &EnvState::get_pressure, &EnvState::set_pressure)
     ;
 
     py::class_<Scenario>(m,
@@ -271,6 +273,16 @@ PYBIND11_MODULE(_PyPartMC, m) {
         "Convert diameter (m) to radius (m)."
     );
 
+    m.def(
+        "loss_rate_dry_dep", &loss_rate_dry_dep, py::return_value_policy::copy,
+        "Compute and return the dry deposition rate for a given particle."
+    );
+
+    m.def(
+        "loss_rate", &loss_rate, py::return_value_policy::copy,
+        "Evaluate a loss rate function."
+    );
+
     m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO);
 
     m.attr("__all__") = py::make_tuple(
@@ -293,6 +305,8 @@ PYBIND11_MODULE(_PyPartMC, m) {
         "sphere_vol2rad",
         "rad2diam",
         "sphere_rad2vol",
-        "diam2rad"
+        "diam2rad",
+        "loss_rate_dry_dep",
+        "loss_rate"
     );
 }

src/scenario.F90

@@ -39,4 +39,58 @@ subroutine f_scenario_from_json(gas_ptr_c, aer_ptr_c, ptr_c) bind(C)
     call spec_file_read_scenario(file, gas_ptr_f, aer_ptr_f, ptr_f)
   end subroutine
 
+  subroutine f_scenario_loss_rate( &
+    scenario_ptr_c, vol, &
+    density, aero_data_ptr_c, &
+    env_state_ptr_c, &
+    rate &
+  ) bind (C)
+
+    type(scenario_t), pointer :: scenario_ptr_f => null()
+    type(aero_data_t), pointer :: aero_data_ptr_f => null()
+    type(env_state_t), pointer :: env_state_ptr_f => null()
+    real(c_double), intent(in) :: vol, density
+    type(c_ptr), intent(in) :: scenario_ptr_c, aero_data_ptr_c, env_state_ptr_c
+    real(c_double), intent(out) :: rate
+
+    call c_f_pointer(scenario_ptr_c, scenario_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)
+
+    rate = scenario_loss_rate( &
+      scenario_ptr_f, &
+      vol, &
+      density, &
+      aero_data_ptr_f, &
+      env_state_ptr_f &
+    )
+
+  end subroutine
+
+  subroutine f_scenario_loss_rate_dry_dep( &
+    vol, &
+    density, &
+    aero_data_ptr_c, &
+    env_state_ptr_c, &
+    rate &
+  ) bind(C)
+
+    type(aero_data_t), pointer :: aero_data_ptr_f => null()
+    type(env_state_t), pointer :: env_state_ptr_f => null()
+    real(c_double), intent(in) :: vol, density
+    type(c_ptr), intent(in) :: aero_data_ptr_c, env_state_ptr_c
+    real(c_double), intent(out) :: rate
+
+    call c_f_pointer(aero_data_ptr_c, aero_data_ptr_f)
+    call c_f_pointer(env_state_ptr_c, env_state_ptr_f)
+
+    rate = scenario_loss_rate_dry_dep( &
+      vol, &
+      density, &
+      aero_data_ptr_f, &
+      env_state_ptr_f &
+    )
+
+  end subroutine
+
 end module

src/scenario.hpp

@@ -9,6 +9,7 @@
 #include "gimmicks.hpp"
 #include "pmc_resource.hpp"
 #include "aero_data.hpp"
+#include "env_state.hpp"
 #include "gas_data.hpp"
 
 extern "C" void f_scenario_ctor(void *ptr) noexcept;
@@ -22,6 +23,21 @@ extern "C" void f_scenario_from_json(
     noexcept
 #endif
 ;
+extern "C" void f_scenario_loss_rate(
+    const void *scenario,
+    const double *vol,
+    const double *density,
+    const void *aero_data,
+    const void *env_state,
+    double *rate
+) noexcept;
+extern "C" void f_scenario_loss_rate_dry_dep(
+    const double *vol,
+    const double *density,
+    const void *aero_data,
+    const void *env_state,
+    double *rate
+) noexcept;
 
 struct Scenario {
     PMCResource ptr;
@@ -49,3 +65,38 @@ struct Scenario {
     }   
 };
 
+double loss_rate(
+    const Scenario &scenario,
+    const double vol,
+    const double density,
+    const AeroData &aero_data,
+    const EnvState &env_state
+) {
+    double rate;
+    f_scenario_loss_rate(
+        &scenario.ptr,
+        &vol,
+        &density,
+        &aero_data.ptr,
+        &env_state.ptr,
+        &rate
+    );
+    return rate;
+}
+
+double loss_rate_dry_dep(
+    const double vol,
+    const double density,
+    const AeroData &aero_data,
+    const EnvState &env_state
+) {
+    double rate;
+    f_scenario_loss_rate_dry_dep(
+        &vol,
+        &density,
+        &aero_data.ptr,
+        &env_state.ptr,
+        &rate
+    );
+    return rate;
+}

tests/test_env_state.py

@@ -44,3 +44,27 @@ def test_dtor():
         # assert
         pass
 
+    @staticmethod
+    def test_height():
+        # arrange
+        sut = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL)
+        value = 1
+
+        # act
+        sut.height = value
+
+        # assert
+        assert value == sut.height
+
+    @staticmethod
+    def test_pressure():
+        # arrange
+        sut = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL)
+        value = 101325
+
+        # act
+        sut.pressure = value
+
+        # assert
+        assert value == sut.pressure
+

tests/test_scenario.py

@@ -4,13 +4,17 @@
 # Authors: https://github.com/open-atmos/PyPartMC/graphs/contributors                              #
 ####################################################################################################
 
+from cmath import nan
 import gc
 import json
 import pytest
+import numpy as np
 import PyPartMC as ppmc
+from PyPartMC import si
 
 from .test_gas_data import GAS_DATA_CTOR_ARG_MINIMAL
 from .test_aero_data import AERO_DATA_CTOR_ARG_MINIMAL
+from .test_env_state import ENV_STATE_CTOR_ARG_MINIMAL
 
 
 SCENARIO_CTOR_ARG_MINIMAL = {
@@ -117,3 +121,59 @@ def test_str():
         # assert
         assert json_actual == SCENARIO_CTOR_ARG_MINIMAL
 
+    @staticmethod
+    @pytest.mark.parametrize("loss_function_params", (
+        {'loss_function': 'none'},
+        {'loss_function': 'constant'},
+        {'loss_function': 'volume'},
+        {'loss_function': 'drydep'},
+        {
+            'loss_function': 'chamber',
+            'chamber_vol': 84.3 * si.m**3,
+            'area_diffuse': 103 * si.m**2,
+            'area_sedi': 12.6 * si.m**2,
+            'prefactor_BL': 0.005 * si.m,
+            'exponent_BL': 0.274
+        }
+    ))
+    def test_loss_rate(loss_function_params:str):
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        env_state = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL)
+        gas_data  = ppmc.GasData(GAS_DATA_CTOR_ARG_MINIMAL)
+        scenario_ctor_arg = {**SCENARIO_CTOR_ARG_MINIMAL, **loss_function_params}
+        scenario = ppmc.Scenario(gas_data, aero_data, scenario_ctor_arg)
+        vol =  (4/3)*np.pi*(1e-6)**3
+        density = 1
+        env_state.height = 1
+        env_state.set_temperature(300)
+        env_state.pressure = 101325
+        aero_data.frac_dim = 3
+        aero_data.prime_radius = 1e-8
+        aero_data.vol_fill_factor = 1
+
+        # act
+        rate = ppmc.loss_rate(scenario, vol, density, aero_data, env_state)
+
+        # assert
+        assert rate is not nan
+
+    @staticmethod
+    def test_loss_rate_dry_dep():
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        env_state = ppmc.EnvState(ENV_STATE_CTOR_ARG_MINIMAL)
+        vol =  (4/3)*np.pi*(1e-6)**3
+        density = 1
+        env_state.height = 1
+        env_state.set_temperature(300)
+        env_state.pressure = 101325
+        aero_data.frac_dim = 3
+        aero_data.prime_radius = 1e-8
+        aero_data.vol_fill_factor = 1
+
+        # act
+        rate = ppmc.loss_rate_dry_dep(vol, density, aero_data, env_state)
+
+        # assert
+        assert rate is not nan