new example based on Grabowski & Pawlowska 2023 (GRL) (#1093)

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

.zenodo.json

@@ -62,6 +62,10 @@
       "affiliation": "California Institute of Technology, Pasadena, CA, USA",
       "name": "Mints, Mikhail"
     },
+    {
+      "affiliation": "University of Warsaw, Poland",
+      "name": "Makulska, Agnieszka"
+    },
     {
       "affiliation": "Jagiellonian University, Kraków, Poland",
       "name": "Olesik, Michael",

PySDM/attributes/impl/mapper.py

@@ -32,6 +32,9 @@
     DryVolumeOrganic,
     OrganicFraction,
 )
+from PySDM.attributes.physics.equilibrium_supersaturation import (
+    EquilibriumSupersaturation,
+)
 from PySDM.attributes.physics.hygroscopicity import Kappa, KappaTimesDryVolume
 from PySDM.attributes.physics.relative_fall_velocity import RelativeFallMomentum
 from PySDM.dynamics.impl.chemistry_utils import AQUEOUS_COMPOUNDS
@@ -93,6 +96,7 @@
     "freezing temperature": lambda _, __: FreezingTemperature,
     "immersed surface area": lambda _, __: ImmersedSurfaceArea,
     "critical supersaturation": lambda _, __: CriticalSupersaturation,
+    "equilibrium supersaturation": lambda _, __: EquilibriumSupersaturation,
     "wet to critical volume ratio": lambda _, __: WetToCriticalVolumeRatio,
 }
 

PySDM/attributes/physics/__init__.py

@@ -6,6 +6,7 @@
 from .critical_volume import CriticalVolume, WetToCriticalVolumeRatio
 from .dry_radius import DryRadius
 from .dry_volume import DryVolume
+from .equilibrium_supersaturation import EquilibriumSupersaturation
 from .heat import Heat
 from .multiplicities import Multiplicities
 from .radius import Radius

PySDM/attributes/physics/equilibrium_supersaturation.py

@@ -0,0 +1,39 @@
+"""
+kappa-Koehler equilibrium supersaturation calculated for actual environment temperature
+"""
+from PySDM.attributes.impl.derived_attribute import DerivedAttribute
+
+
+class EquilibriumSupersaturation(DerivedAttribute):
+    def __init__(self, builder):
+        self.r_wet = builder.get_attribute("radius")
+        self.v_wet = builder.get_attribute("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="equilibrium supersaturation",
+            dependencies=(self.kappa, self.v_dry, self.f_org, self.r_wet),
+        )
+
+    def recalculate(self):
+        if len(self.particulator.environment["T"]) != 1:
+            raise NotImplementedError()
+        temperature = self.particulator.environment["T"][0]
+        rd3 = self.v_dry.data.data / self.formulae.constants.PI_4_3
+        sgm = self.formulae.surface_tension.sigma(
+            temperature,
+            self.v_wet.data.data,
+            self.v_dry.data.data,
+            self.f_org.data.data,
+        )
+
+        self.data.data[:] = self.formulae.hygroscopicity.RH_eq(
+            self.r_wet.data.data,
+            T=temperature,
+            kp=self.kappa.data.data,
+            rd3=rd3,
+            sgm=sgm,
+        )

PySDM/products/impl/activation_filtered_product.py

@@ -0,0 +1,30 @@
+"""
+common base class for products filtering droplets based on their activation state
+"""
+import numpy as np
+
+
+class _ActivationFilteredProduct:
+    def __init__(
+        self,
+        *,
+        count_unactivated: bool,
+        count_activated: bool,
+    ):
+        self.__filter_attr = "wet to critical volume ratio"
+        self.__filter_range = [0, np.inf]
+        if not count_activated:
+            self.__filter_range[1] = 1
+        if not count_unactivated:
+            self.__filter_range[0] = 1
+
+    def impl(self, *, attr, rank):
+        getattr(self, "_download_moment_to_buffer")(
+            attr=attr,
+            rank=rank,
+            filter_attr=self.__filter_attr,
+            filter_range=self.__filter_range,
+        )
+
+    def register(self, builder):
+        builder.request_attribute(self.__filter_attr)

PySDM/products/impl/concentration_product.py

@@ -8,6 +8,10 @@
 
 class ConcentrationProduct(MomentProduct):
     def __init__(self, *, unit: str, name: str, specific: bool, stp: bool):
+        """
+        `stp` toggles expressing the concentration in terms of standard temperature
+        and pressure conditions (ground level of the ICAO standard atmosphere, zero humidity)
+        """
         super().__init__(unit=unit, name=name)
         self.specific = specific
         self.stp = stp

PySDM/products/size_spectral/__init__.py

@@ -7,8 +7,14 @@
 )
 from .effective_radius import EffectiveRadius
 from .mean_radius import MeanRadius
+from .mean_radius_activated import ActivatedMeanRadius
+from .mean_volume_radius import MeanVolumeRadius
 from .number_size_spectrum import NumberSizeSpectrum
 from .particle_concentration import ParticleConcentration, ParticleSpecificConcentration
+from .particle_concentration_activated import (
+    ActivatedParticleConcentration,
+    ActivatedParticleSpecificConcentration,
+)
 from .particle_size_spectrum import (
     ParticleSizeSpectrumPerMass,
     ParticleSizeSpectrumPerVolume,
@@ -19,6 +25,11 @@
 from .radius_binned_number_averaged_terminal_velocity import (
     RadiusBinnedNumberAveragedTerminalVelocity,
 )
+from .size_standard_deviation import (
+    AreaStandardDeviation,
+    RadiusStandardDeviation,
+    VolumeStandardDeviation,
+)
 from .total_particle_concentration import TotalParticleConcentration
 from .total_particle_specific_concentration import TotalParticleSpecificConcentration
 from .water_mixing_ratio import WaterMixingRatio

PySDM/products/size_spectral/mean_radius.py

@@ -1,14 +1,29 @@
 """
 mean radius of particles within a grid cell (optionally restricted to a given size range)
 """
+import numpy as np
+
 from PySDM.products.impl.moment_product import MomentProduct
 
 
 class MeanRadius(MomentProduct):
-    def __init__(self, name=None, unit="m"):
+    def __init__(
+        self,
+        name=None,
+        unit="m",
+        radius_range=(0, np.inf),
+    ):
+        self.radius_range = radius_range
         super().__init__(name=name, unit=unit)
 
     def _impl(self, **kwargs):
-        self._download_moment_to_buffer(attr="volume", rank=1 / 3)
+        self._download_moment_to_buffer(
+            attr="volume",
+            rank=1 / 3,
+            filter_range=(
+                self.formulae.trivia.volume(self.radius_range[0]),
+                self.formulae.trivia.volume(self.radius_range[1]),
+            ),
+        )
         self.buffer[:] /= self.formulae.constants.PI_4_3 ** (1 / 3)
         return self.buffer

PySDM/products/size_spectral/mean_radius_activated.py

@@ -0,0 +1,24 @@
+"""
+mean radius of particles within a grid cell, for activated, unactivated or both
+"""
+from PySDM.products.impl.activation_filtered_product import _ActivationFilteredProduct
+from PySDM.products.impl.moment_product import MomentProduct
+
+
+class ActivatedMeanRadius(MomentProduct, _ActivationFilteredProduct):
+    def __init__(
+        self, count_unactivated: bool, count_activated: bool, name=None, unit="m"
+    ):
+        MomentProduct.__init__(self, name=name, unit=unit)
+        _ActivationFilteredProduct.__init__(
+            self, count_activated=count_activated, count_unactivated=count_unactivated
+        )
+
+    def register(self, builder):
+        for base_class in (_ActivationFilteredProduct, MomentProduct):
+            base_class.register(self, builder)
+
+    def _impl(self, **kwargs):
+        _ActivationFilteredProduct.impl(self, attr="volume", rank=1 / 3)
+        self.buffer[:] /= self.formulae.constants.PI_4_3 ** (1 / 3)
+        return self.buffer

PySDM/products/size_spectral/mean_volume_radius.py

@@ -0,0 +1,24 @@
+"""
+mean volume radius of particles within a grid cell, for activated, unactivated or both
+"""
+
+from PySDM.products.impl.activation_filtered_product import _ActivationFilteredProduct
+from PySDM.products.impl.moment_product import MomentProduct
+
+
+class MeanVolumeRadius(MomentProduct, _ActivationFilteredProduct):
+    def __init__(
+        self, count_unactivated: bool, count_activated: bool, name=None, unit="m"
+    ):
+        MomentProduct.__init__(self, name=name, unit=unit)
+        _ActivationFilteredProduct.__init__(
+            self, count_activated=count_activated, count_unactivated=count_unactivated
+        )
+
+    def register(self, builder):
+        for base_class in (_ActivationFilteredProduct, MomentProduct):
+            base_class.register(self, builder)
+
+    def _impl(self, **kwargs):
+        _ActivationFilteredProduct.impl(self, attr="volume", rank=1)
+        return self.formulae.trivia.radius(self.buffer[:])

PySDM/products/size_spectral/particle_concentration_activated.py

@@ -0,0 +1,46 @@
+"""
+concentration of particles within a grid cell (either per-volume of per-mass-of-dry air),
+ activated, unactivated or both
+"""
+
+from PySDM.products.impl.activation_filtered_product import _ActivationFilteredProduct
+from PySDM.products.impl.concentration_product import ConcentrationProduct
+
+
+class ActivatedParticleConcentration(ConcentrationProduct, _ActivationFilteredProduct):
+    # pylint: disable=too-many-arguments
+    def __init__(
+        self,
+        *,
+        count_unactivated: bool,
+        count_activated: bool,
+        specific=False,
+        stp=False,
+        name=None,
+        unit="m^-3",
+    ):
+        ConcentrationProduct.__init__(
+            self, name=name, unit=unit, specific=specific, stp=stp
+        )
+        _ActivationFilteredProduct.__init__(
+            self, count_activated=count_activated, count_unactivated=count_unactivated
+        )
+
+    def register(self, builder):
+        for base_class in (_ActivationFilteredProduct, ConcentrationProduct):
+            base_class.register(self, builder)
+
+    def _impl(self, **kwargs):
+        _ActivationFilteredProduct.impl(self, attr="volume", rank=0)
+        return ConcentrationProduct._impl(self, **kwargs)
+
+
+class ActivatedParticleSpecificConcentration(ActivatedParticleConcentration):
+    def __init__(self, count_unactivated, count_activated, name=None, unit="kg^-1"):
+        super().__init__(
+            count_unactivated=count_unactivated,
+            count_activated=count_activated,
+            specific=True,
+            name=name,
+            unit=unit,
+        )

PySDM/products/size_spectral/size_standard_deviation.py

@@ -0,0 +1,68 @@
+"""
+standard deviation of radius/area/volume of particles within a grid cell,
+for activated, unactivated or both
+"""
+import numpy as np
+
+from PySDM.products.impl.activation_filtered_product import _ActivationFilteredProduct
+from PySDM.products.impl.moment_product import MomentProduct
+
+
+class _SizeStandardDeviation(MomentProduct, _ActivationFilteredProduct):
+    # pylint: disable=too-many-arguments
+    def __init__(
+        self,
+        count_unactivated: bool,
+        count_activated: bool,
+        name=None,
+        unit="m",
+        attr="radius",
+    ):
+        self.attr = attr
+        MomentProduct.__init__(self, name=name, unit=unit)
+        _ActivationFilteredProduct.__init__(
+            self, count_activated=count_activated, count_unactivated=count_unactivated
+        )
+
+    def register(self, builder):
+        builder.request_attribute(self.attr)
+        for base_class in (_ActivationFilteredProduct, MomentProduct):
+            base_class.register(self, builder)
+
+    def _impl(self, **kwargs):
+        _ActivationFilteredProduct.impl(self, attr=self.attr, rank=1)
+        tmp = np.empty_like(self.buffer)
+        tmp[:] = -self.buffer**2
+        _ActivationFilteredProduct.impl(self, attr=self.attr, rank=2)
+        tmp[:] += self.buffer
+        tmp[:] = np.sqrt(tmp)
+        return tmp
+
+
+RadiusStandardDeviation = _SizeStandardDeviation
+
+
+class AreaStandardDeviation(_SizeStandardDeviation):
+    def __init__(
+        self, *, name=None, unit="m^2", count_activated: bool, count_unactivated: bool
+    ):
+        super().__init__(
+            name=name,
+            unit=unit,
+            count_activated=count_activated,
+            count_unactivated=count_unactivated,
+            attr="area",
+        )
+
+
+class VolumeStandardDeviation(_SizeStandardDeviation):
+    def __init__(
+        self, *, name=None, unit="m^3", count_activated: bool, count_unactivated: bool
+    ):
+        super().__init__(
+            name=name,
+            unit=unit,
+            count_activated=count_activated,
+            count_unactivated=count_unactivated,
+            attr="volume",
+        )

examples/PySDM_examples/Grabowski_and_Pawlowska_2023/__init__.py

@@ -0,0 +1,3 @@
+# pylint: disable=invalid-name
+from .settings import Settings
+from .simulation import Simulation