Wall loss example

docs/_toc.yml

@@ -23,7 +23,8 @@ parts:
     - file: examples/streamlake/stream_stats_size_distribution_part2
   - file: examples/wall_loss_section
     sections:
-    - file: examples/chamber_smps_data
+    - file: examples/wall_loss/chamber_smps_data
+    - file: examples/wall_loss/wall_loss_forward_simulation
   - file: examples/lagrangian/lagrangian_intro
     sections:
     - file: examples/lagrangian/basic_lagrangian_box

particula/dynamics.py

@@ -1,6 +1,6 @@
 """ statics -> dynamics
 """
-from scipy.integrate import odeint
+from scipy.integrate import odeint, solve_ivp
 
 from particula.rates import Rates
 
@@ -12,11 +12,13 @@ class Solver(Rates):
     def __init__(
         self,
         time_span=None,
-        do_coag=True,
-        do_cond=True,
-        do_nucl=True,
+        do_coagulation=True,
+        do_condensation=True,
+        do_nucleation=True,
+        do_dilution=False,
+        do_wall_loss=False,
         **kwargs
-    ):
+    ):  # pylint: disable=too-many-arguments
         """ constructor
         """
         super().__init__(**kwargs)
@@ -25,9 +27,11 @@ def __init__(
             raise ValueError("You must provide a time span!")
 
         self.time_span = time_span
-        self.do_coag = do_coag
-        self.do_cond = do_cond
-        self.do_nucl = do_nucl
+        self.do_coagulation = do_coagulation
+        self.do_condensation = do_condensation
+        self.do_nucleation = do_nucleation
+        self.do_dilution = do_dilution
+        self.do_wall_loss = do_wall_loss
 
     def _ode_func(self, _nums, _,):
         """ ode_func
@@ -39,17 +43,37 @@ def _ode_func(self, _nums, _,):
         )
 
         return self.sum_rates(
-            coagulation=self.do_coag,
-            condensation=self.do_cond,
-            nucleation=self.do_nucl,
+            coagulation=self.do_coagulation,
+            condensation=self.do_condensation,
+            nucleation=self.do_nucleation,
+            dilution=self.do_dilution,
+            wall_loss=self.do_wall_loss,
         ).m
 
-    def solution(self):
+    def solution(
+            self,
+            method='odeint',
+            **kwargs_ode
+    ):
         """ solve the equation
         """
+        if method == 'odeint':
+            return odeint(
+                func=self._ode_func,
+                y0=self.particle.particle_distribution().m,
+                t=self.time_span,
+                **kwargs_ode
+            )*self.particle.particle_distribution().u
+
+        if method == 'solve_ivp':
+            if 'method' not in kwargs_ode:
+                kwargs_ode['method'] = 'BDF'  # Choose an appropriate method
+            return solve_ivp(
+                fun=self._ode_func,
+                t_span=(self.time_span[0], self.time_span[-1]),
+                y0=self.particle.particle_distribution().m,
+                t_eval=self.time_span,
+                **kwargs_ode
+            ).y.T*self.particle.particle_distribution().u
 
-        return odeint(
-            func=self._ode_func,
-            y0=self.particle.particle_distribution().m,
-            t=self.time_span,
-        )*self.particle.particle_distribution().u
+        raise ValueError("Invalid method!")

particula/particle.py

@@ -23,6 +23,9 @@
 from particula.util.rms_speed import cbar
 from particula.util.slip_correction import scf
 from particula.util.vapor_flux import phi
+from particula.util.settling_velocity import psv
+from particula.util.aerodynamic_mobility import pam
+from particula.util.diffusion_coefficient import pdc
 from particula.util.wall_loss import wlc
 from particula.util.kelvin_correction import kelvin_term
 from particula.vapor import Vapor
@@ -221,6 +224,33 @@ def friction_factor(self):
             scf=self.slip_correction_factor(),
         )
 
+    def settling_velocity(self):
+        """ Returns a particle's settling velocity.
+        """
+        return psv(
+            rad=self.particle_radius,
+            den=self.particle_density,
+            scf_val=self.slip_correction_factor(),
+            vis_val=self.dynamic_viscosity(),
+        )
+
+    def aerodynamic_mobility(self):
+        """ Returns a particle's aerodynamic mobility.
+        """
+        return pam(
+            radius=self.particle_radius,
+            scf_val=self.slip_correction_factor(),
+            vis_val=self.dynamic_viscosity(),
+        )
+
+    def diffusion_coefficient(self):
+        """ Returns a particle's diffusion coefficient.
+        """
+        return pdc(
+            temperature=self.temperature,
+            pam_val=self.aerodynamic_mobility(),
+        )
+
 
 class ParticleCondensation(ParticleInstances):
     """ calculate some condensation stuff
@@ -332,17 +362,25 @@ def __init__(self, **kwargs):
         self.wall_loss_approximation = str(
             kwargs.get("wall_loss_approximation", "none")
         )
-        self.chamber_radius = in_length(
-            kwargs.get("chamber_radius", 3)
+        self.chamber_dimension = in_length(
+            kwargs.get("chamber_dimension", 1)
+        )
+        self.chamber_ktp_value = in_handling(
+            kwargs.get("chamber_ktp_value", 0.1),
+            u.s**-1
         )
         self.kwargs = kwargs
 
     def wall_loss_coefficient(self):
         """ Returns a particle's wall loss coefficient.
         """
         return wlc(
-            approx=self.wall_loss_approximation
-        )
+            approx=self.wall_loss_approximation,
+            ktp_value=self.chamber_ktp_value,
+            diffusion_coefficient_value=self.diffusion_coefficient(),
+            dimension=self.chamber_dimension,
+            settling_velocity_value=self.settling_velocity(),
+            )
 
 
 class Particle(ParticleWallLoss):

particula/rates.py

@@ -108,33 +108,47 @@ def wall_loss_rate(self):
             self.particle_distribution
         )
 
-    def sum_rates(self, coagulation=True, condensation=True, nucleation=True):
+    def sum_rates(
+            self,
+            coagulation=True,
+            condensation=True,
+            nucleation=True,
+            dilution=False,
+            wall_loss=False,
+    ):  # pylint: disable=too-many-arguments
         """Sum rates, with options to disable individual rate terms.
 
-        Parameters
+        Args:
         ----------
         coagulation : bool, optional
             does the coagulation calcuation, by default True
         condensation : bool, optional
             does the condensation calculation, by default True
         nucleation : bool, optional
             does the nucleation calculation, by default True
+        dilution : bool, optional
+            does the dilution calculation, by default False
+        wall_loss : bool, optional
+            does the wall loss calculation, by default False
 
-        TODO: add wall and dilution loss rates
         """
 
         # Define a dictionary that maps option names to rate functions
         rates = {
             'coagulation': self.coagulation_rate,
             'condensation': self.condensation_growth_rate,
-            'nucleation': self.nucleation_rate
+            'nucleation': self.nucleation_rate,
+            'dilution': self.dilution_rate,
+            'wall_loss': self.wall_loss_rate,
         }
 
         # Define a dictionary that maps option names to their Boolean values
         options = {
             'coagulation': coagulation,
             'condensation': condensation,
-            'nucleation': nucleation
+            'nucleation': nucleation,
+            'dilution': dilution,
+            'wall_loss': wall_loss,
         }
         # Return the sum of the rates that are enabled
         return sum(

particula/util/aerodynamic_mobility.py

@@ -18,6 +18,8 @@ def pam(
     rad = in_radius(radius)
     scf_val = in_scalar(
         scf_val) if scf_val is not None else scf(radius=rad, **kwargs)
+
+    # dyn_vis has no radius, is this a bug waiting to happen?
     vis_val = in_viscosity(
         vis_val) if vis_val is not None else dyn_vis(radius=rad, **kwargs)
 

particula/util/tests/wall_loss_test.py

@@ -0,0 +1,53 @@
+"""
+Test the wall loss functions.
+"""
+
+import numpy as np
+
+from particula.util.wall_loss import (
+    spherical_wall_loss_coefficient,
+    rectangle_wall_loss)
+
+
+def test_spherical_wall_loss_coefficient():
+    """test the spherical wall loss coefficient function"""
+    # Test case with hypothetical values
+    ktp_value = 1.5  # Example value
+    diffusion_coefficient_value = 0.01  # Example value
+    settling_velocity_value = 0.05  # Example value
+    chamber_radius = 2.0  # Example value
+
+    expected_output = 0.11816  # This is a hypothetical value.
+
+    # Call the function with the test case
+    calculated_output = spherical_wall_loss_coefficient(
+        ktp_value,
+        diffusion_coefficient_value,
+        settling_velocity_value,
+        chamber_radius
+    )
+
+    # Assertion to check if the calculated output matches the expected output
+    assert np.isclose(calculated_output, expected_output, rtol=1e-4)
+
+
+def test_rectangle_wall_loss():
+    """test the rectangular wall loss coefficient function"""
+    # Test case with hypothetical values
+    ktp_value = 1.5  # Example value
+    diffusion_coefficient_value = 0.01  # Example value
+    settling_velocity_value = 0.05  # Example value
+    dimension = [1, 1, 1]  # Example value
+
+    expected_output = 0.47312  # This is a hypothetical value.
+
+    # Call the function with the test case
+    calculated_output = rectangle_wall_loss(
+        ktp_value,
+        diffusion_coefficient_value,
+        settling_velocity_value,
+        dimension
+    )
+
+    # Assertion to check if the calculated output matches the expected output
+    assert np.isclose(calculated_output, expected_output, rtol=1e-4)