Merge pull request #657 from slayoo/dev

adding critical supersaturation attribute and activable fraction product

PySDM/attributes/impl/mapper.py

@@ -7,6 +7,7 @@
 from PySDM.attributes.ice import FreezingTemperature, ImmersedSurfaceArea
 from PySDM.attributes.numerics import CellID, CellOrigin, PositionInCell
 from PySDM.attributes.chemistry import MoleAmount, Concentration, pH, HydrogenIonConcentration
+from PySDM.attributes.physics.critical_supersaturation import CriticalSupersaturation
 from PySDM.physics.aqueous_chemistry.support import AQUEOUS_COMPOUNDS
 from PySDM.physics.surface_tension import Constant
 
@@ -15,14 +16,20 @@
     'volume': lambda _: Volume,
     'dry volume organic': lambda dynamics: (
         DummyAttributeImpl('dry volume organic')
-        if isinstance(dynamics['Condensation'].particulator.formulae.surface_tension, Constant)
+        if 'Condensation' in dynamics and isinstance(
+            dynamics['Condensation'].particulator.formulae.surface_tension,
+            Constant
+        )
         else DryVolumeOrganic
     ),
     'dry volume': lambda dynamics:
     DryVolumeDynamic if 'AqueousChemistry' in dynamics else DryVolume,
     'dry volume organic fraction': lambda dynamics: (
         DummyAttributeImpl('dry volume organic fraction')
-        if isinstance(dynamics['Condensation'].particulator.formulae.surface_tension, Constant)
+        if 'Condensation' in dynamics and isinstance(
+            dynamics['Condensation'].particulator.formulae.surface_tension,
+            Constant
+        )
         else OrganicFraction
     ),
     'kappa times dry volume': lambda _: KappaTimesDryVolume,
@@ -43,7 +50,8 @@
     'pH': lambda _: pH,
     'conc_H': lambda _: HydrogenIonConcentration,
     'freezing temperature': lambda _: FreezingTemperature,
-    'immersed surface area': lambda _: ImmersedSurfaceArea
+    'immersed surface area': lambda _: ImmersedSurfaceArea,
+    'critical supersaturation': lambda _: CriticalSupersaturation
 }
 
 

PySDM/attributes/physics/__init__.py

@@ -7,3 +7,4 @@
 from .temperature import Temperature
 from .heat import Heat
 from .critical_volume import CriticalVolume
+from .critical_supersaturation import CriticalSupersaturation

PySDM/attributes/physics/critical_supersaturation.py

@@ -0,0 +1,36 @@
+from PySDM.attributes.impl.derived_attribute import DerivedAttribute
+from PySDM.physics import constants as const
+
+
+class CriticalSupersaturation(DerivedAttribute):
+    def __init__(self, builder):
+        self.v_crit = builder.get_attribute('critical volume')
+        self.v_dry = builder.get_attribute('dry volume')
+        self.kappa = builder.get_attribute('kappa')
+        self.f_org = builder.get_attribute('dry volume organic fraction')
+
+        super().__init__(
+            builder=builder,
+            name='critical supersaturation',
+            dependencies=(self.v_crit, self.kappa, self.v_dry, self.f_org)
+        )
+        self.formulae = builder.particulator.formulae
+        self.environment = builder.particulator.environment
+
+    def recalculate(self):
+        if len(self.environment['T']) != 1:
+            raise NotImplementedError()
+        T = self.environment['T'][0]
+        r_cr = self.formulae.trivia.radius(self.v_crit.data.data)
+        rd3 = self.v_dry.data.data / const.pi_4_3
+        sgm = self.formulae.surface_tension.sigma(
+            T, self.v_crit.data.data, self.v_dry.data.data, self.f_org.data.data
+        )
+
+        self.data.data[:] = self.formulae.hygroscopicity.RH_eq(
+            r_cr,
+            T=T,
+            kp=self.kappa.data.data,
+            rd3=rd3,
+            sgm=sgm
+        )

PySDM/attributes/physics/dry_radius.py

@@ -9,6 +9,5 @@ def __init__(self, builder):
         super().__init__(builder, name='dry radius', dependencies=dependencies)
 
     def recalculate(self):
-        self.data.idx = self.volume_dry.data.idx
         self.data.product(self.volume_dry.get(), 1/const.pi_4_3)
         self.data **= 1/3

PySDM/products/dynamics/condensation/__init__.py

@@ -1,3 +1,4 @@
 from .condensation_timestep import CondensationTimestepMin, CondensationTimestepMax
 from .event_rates import RipeningRate, ActivatingRate, DeactivatingRate
 from .peak_supersaturation import PeakSupersaturation
+from .activable_fraction import ActivableFraction

PySDM/products/dynamics/condensation/activable_fraction.py

@@ -0,0 +1,29 @@
+from ...product import MomentProduct
+
+
+class ActivableFraction(MomentProduct):
+    def __init__(self):
+        super().__init__(
+            name="activable fraction",
+            unit="1",
+            description=""
+        )
+
+    def register(self, builder):
+        super().register(builder)
+        builder.request_attribute('critical supersaturation')
+
+    def get(self, S_max):
+        self.download_moment_to_buffer(
+            'volume',
+            rank=0,
+            filter_range=(0, 1 + S_max / 100),
+            filter_attr='critical supersaturation'
+        )
+        frac = self.buffer.copy()
+        self.download_moment_to_buffer(
+            'volume',
+            rank=0
+        )
+        frac /= self.buffer
+        return frac

PySDM/products/dynamics/freezing/ice_water_content.py

@@ -1,6 +1,7 @@
+import numpy as np
 from ...product import MomentProduct
 from ....physics import constants as const
-import numpy as np
+
 
 class IceWaterContent(MomentProduct):
 
@@ -24,4 +25,4 @@ def get(self):
             self.download_to_buffer(self.particulator.environment['rhod'])
             result[:] /= self.buffer
             const.convert_to(result, const.si.gram / const.si.kilogram)
-        return result
+        return result

tests/unit_tests/attributes/physics/test_critical_supersaturation.py

@@ -0,0 +1,35 @@
+import numpy as np
+from PySDM.products.dynamics.condensation import ActivableFraction
+from PySDM import Builder
+from PySDM.backends import CPU
+from PySDM.environments import Box
+from PySDM.physics import si
+
+
+def test_critical_supersaturation():
+    # arrange
+    T = 300 * si.K
+    n_sd = 100
+    S_max = .01
+    vdry = np.linspace(.001, 1, n_sd) * si.um**3
+
+    builder = Builder(n_sd=n_sd, backend=CPU())
+    env = Box(dt=np.nan, dv=np.nan)
+    builder.set_environment(env)
+    env['T'] = T
+    particulator = builder.build(
+        attributes={
+            'n': np.ones(n_sd),
+            'volume': np.linspace(.01, 10, n_sd) * si.um**3,
+            'dry volume': vdry,
+            'kappa times dry volume': .9 * vdry,
+            'dry volume organic': np.zeros(n_sd)
+        },
+        products=[ActivableFraction()]
+    )
+
+    # act
+    AF = particulator.products['activable fraction'].get(S_max)
+
+    # assert
+    assert 0 < AF < 1