homogeneous nucleation of liquid droplets & expansion chamber example

.github/workflows/tests.yml

@@ -84,7 +84,11 @@ jobs:
       matrix:
         platform: [ubuntu-24.04, macos-13, macos-14, windows-latest]
         python-version: ["3.9", "3.12"]
+<<<<<<< HEAD
+        test-suite: ["unit_tests/!(dynamics)", "unit_tests/dynamics/!(collisions)", "unit_tests/dynamics/collisions", "smoke_tests/no_env", "smoke_tests/box", "smoke_tests/parcel_a", "smoke_tests/parcel_b", "smoke_tests/parcel_c", "smoke_tests/parcel_d", "smoke_tests/kinematic_1d", "smoke_tests/kinematic_2d", "smoke_tests/chamber", "tutorials_tests"]
+=======
         test-suite: ["unit_tests/!(dynamics)", "unit_tests/dynamics/!(condensation|collisions)", "unit_tests/dynamics/condensation", "unit_tests/dynamics/collisions","smoke_tests/no_env", "smoke_tests/box", "smoke_tests/parcel_a", "smoke_tests/parcel_b", "smoke_tests/parcel_c", "smoke_tests/parcel_d", "smoke_tests/kinematic_1d", "smoke_tests/kinematic_2d", "tutorials_tests"]
+>>>>>>> oa/main
         exclude:
           - platform: "macos-14"
             python-version: "3.9"

.gitignore

@@ -4,9 +4,15 @@ PySDM_examples/utils/temporary_files/*
 *.pkl
 *.vtu
 *.vts
+<<<<<<< HEAD
+*.png
+*.csv
+*.xlsx
+=======
 *.svg
 *.gif
 *.ipynb.badges.md
+>>>>>>> oa/main
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

PySDM/backends/impl_numba/methods/__init__.py

@@ -12,5 +12,5 @@
 from .pair_methods import PairMethods
 from .physics_methods import PhysicsMethods
 from .terminal_velocity_methods import TerminalVelocityMethods
-from .seeding_methods import SeedingMethods
+from .spawning_methods import SpawningMethods
 from .deposition_methods import DepositionMethods

PySDM/backends/impl_numba/methods/seeding_methods.py → PySDM/backends/impl_numba/methods/spawning_methods.py

@@ -1,4 +1,4 @@
-"""CPU implementation of backend methods for particle injections"""
+"""CPU implementation of backend methods for particle spawning"""
 
 from functools import cached_property
 
@@ -7,62 +7,62 @@
 from PySDM.backends.impl_common.backend_methods import BackendMethods
 
 
-class SeedingMethods(BackendMethods):  # pylint: disable=too-few-public-methods
+class SpawningMethods(BackendMethods):  # pylint: disable=too-few-public-methods
     @cached_property
-    def _seeding(self):
+    def _spawning(self):
         @numba.njit(**{**self.default_jit_flags, "parallel": False})
         def body(  # pylint: disable=too-many-arguments
             idx,
             multiplicity,
             extensive_attributes,
-            seeded_particle_index,
-            seeded_particle_multiplicity,
-            seeded_particle_extensive_attributes,
-            number_of_super_particles_to_inject: int,
+            spawned_particle_index,
+            spawned_particle_multiplicity,
+            spawned_particle_extensive_attributes,
+            number_of_super_particles_to_spawn: int,
         ):
             number_of_super_particles_already_injected = 0
             # TODO #1387 start enumerating from the end of valid particle set
             for i, mult in enumerate(multiplicity):
                 if (
-                    number_of_super_particles_to_inject
+                    number_of_super_particles_to_spawn
                     == number_of_super_particles_already_injected
                 ):
                     break
                 if mult == 0:
                     idx[i] = -1
-                    s = seeded_particle_index[
+                    s = spawned_particle_index[
                         number_of_super_particles_already_injected
                     ]
                     number_of_super_particles_already_injected += 1
-                    multiplicity[i] = seeded_particle_multiplicity[s]
+                    multiplicity[i] = spawned_particle_multiplicity[s]
                     for a in range(len(extensive_attributes)):
                         extensive_attributes[a, i] = (
-                            seeded_particle_extensive_attributes[a, s]
+                            spawned_particle_extensive_attributes[a, s]
                         )
             assert (
-                number_of_super_particles_to_inject
+                number_of_super_particles_to_spawn
                 == number_of_super_particles_already_injected
             )
 
         return body
 
-    def seeding(
+    def spawning(
         self,
         *,
         idx,
         multiplicity,
         extensive_attributes,
-        seeded_particle_index,
-        seeded_particle_multiplicity,
-        seeded_particle_extensive_attributes,
-        number_of_super_particles_to_inject: int,
+        spawned_particle_index,
+        spawned_particle_multiplicity,
+        spawned_particle_extensive_attributes,
+        number_of_super_particles_to_spawn: int,
     ):
-        self._seeding(
+        self._spawning(
             idx=idx.data,
             multiplicity=multiplicity.data,
             extensive_attributes=extensive_attributes.data,
-            seeded_particle_index=seeded_particle_index.data,
-            seeded_particle_multiplicity=seeded_particle_multiplicity.data,
-            seeded_particle_extensive_attributes=seeded_particle_extensive_attributes.data,
-            number_of_super_particles_to_inject=number_of_super_particles_to_inject,
+            spawned_particle_index=spawned_particle_index.data,
+            spawned_particle_multiplicity=spawned_particle_multiplicity.data,
+            spawned_particle_extensive_attributes=spawned_particle_extensive_attributes.data,
+            number_of_super_particles_to_spawn=number_of_super_particles_to_spawn,
         )

PySDM/backends/numba.py

@@ -28,7 +28,7 @@ class Numba(  # pylint: disable=too-many-ancestors,duplicate-code
     methods.DisplacementMethods,
     methods.TerminalVelocityMethods,
     methods.IsotopeMethods,
-    methods.SeedingMethods,
+    methods.SpawningMethods,
     methods.DepositionMethods,
 ):
     Storage = ImportedStorage
@@ -78,5 +78,5 @@ def __init__(self, formulae=None, double_precision=True, override_jit_flags=None
         methods.DisplacementMethods.__init__(self)
         methods.TerminalVelocityMethods.__init__(self)
         methods.IsotopeMethods.__init__(self)
-        methods.SeedingMethods.__init__(self)
+        methods.SpawningMethods.__init__(self)
         methods.DepositionMethods.__init__(self)

PySDM/dynamics/__init__.py

@@ -16,4 +16,5 @@
 from PySDM.dynamics.freezing import Freezing
 from PySDM.dynamics.relaxed_velocity import RelaxedVelocity
 from PySDM.dynamics.seeding import Seeding
+from PySDM.dynamics.homogeneous_liquid_nucleation import HomogeneousLiquidNucleation
 from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce

PySDM/dynamics/homogeneous_liquid_nucleation.py

@@ -0,0 +1,72 @@
+"""nucleation of droplets in CCN-void air that spawns new super-droplets"""
+
+import numpy as np
+from PySDM.dynamics.impl import register_dynamic
+from PySDM.dynamics.impl import SuperParticleSpawningDynamic
+
+
+@register_dynamic()
+class HomogeneousLiquidNucleation(SuperParticleSpawningDynamic):
+    def __init__(self):
+        self.particulator = None
+        self.formulae = None
+        self.index = None
+
+    def register(self, builder):
+        self.particulator = builder.particulator
+        self.formulae = builder.formulae
+        self.index = self.particulator.Index.identity_index(1)
+
+    def __call__(self):
+        env = {
+            k: self.particulator.environment[k].to_ndarray()[0] for k in ("T", "RH")
+        }  # TODO #1492: >0D
+
+        if env["RH"] < 1:
+            return
+
+        e_s = self.formulae.saturation_vapour_pressure.pvs_water(env["T"])
+        j = self.formulae.homogeneous_liquid_nucleation_rate.j_liq_homo(
+            env["T"], env["RH"], e_s
+        )
+
+        # TODO #1492: take care of cases where round yields zero -> MC sampling?
+        new_sd_multiplicity = round(
+            j * self.particulator.environment.mesh.dv * self.particulator.dt
+        )
+
+        if new_sd_multiplicity > 0:
+            r_wet = self.formulae.homogeneous_liquid_nucleation_rate.r_liq_homo(
+                env["T"], env["RH"]
+            )
+            v_wet = self.formulae.trivia.volume(radius=r_wet)
+            new_sd_extensive_attributes = {
+                "signed water mass": (v_wet * self.formulae.constants.rho_w,),
+                "dry volume": (0,),
+                "kappa times dry volume": (0,),
+            }
+
+            # TODO #1492: to be done once, not every time we spawn
+            self.check_extensive_attribute_keys(
+                particulator_attributes=self.particulator.attributes,
+                spawned_attributes=new_sd_extensive_attributes,
+            )
+
+            # TODO #1492: allocate once, reuse
+            new_sd_extensive_attributes = self.particulator.IndexedStorage.from_ndarray(
+                self.index,
+                np.asarray(list(new_sd_extensive_attributes.values())),
+            )
+            # TODO #1492: ditto
+            new_sd_multiplicity = self.particulator.IndexedStorage.from_ndarray(
+                self.index,
+                np.asarray((new_sd_multiplicity,), dtype=np.int64),
+            )
+
+            self.particulator.spawn(
+                spawned_particle_index=self.index,
+                number_of_super_particles_to_spawn=1,
+                spawned_particle_multiplicity=new_sd_multiplicity,
+                spawned_particle_extensive_attributes=new_sd_extensive_attributes,
+            )
+            # TODO #1492: subtract the water mass from ambient vapour

PySDM/dynamics/impl/__init__.py

@@ -1,3 +1,4 @@
 """stuff not intended to be imported from user code"""
 
 from .register_dynamic import register_dynamic
+from .super_particle_spawning_dynamic import SuperParticleSpawningDynamic

PySDM/dynamics/impl/super_particle_spawning_dynamic.py

@@ -0,0 +1,23 @@
+"""common code for dynamics that spawn new super-droplets"""
+
+from PySDM.impl.particle_attributes import ParticleAttributes
+
+
+class SuperParticleSpawningDynamic:  # pylint: disable=too-few-public-methods
+    """base class for dynamics that spawn new super-droplets"""
+
+    @staticmethod
+    def check_extensive_attribute_keys(
+        particulator_attributes: ParticleAttributes,
+        spawned_attributes: dict,
+    ):
+        """checks if keys (and their order) in user-supplied `spawned_attributes` match
+        attributes used in the `particulator_attributes`"""
+        if tuple(particulator_attributes.get_extensive_attribute_keys()) != tuple(
+            spawned_attributes.keys()
+        ):
+            raise ValueError(
+                f"extensive attributes ({spawned_attributes.keys()})"
+                " do not match those used in particulator"
+                f" ({particulator_attributes.get_extensive_attribute_keys()})"
+            )

PySDM/dynamics/seeding.py

@@ -8,10 +8,11 @@
 
 from PySDM.dynamics.impl import register_dynamic
 from PySDM.initialisation import discretise_multiplicities
+from PySDM.dynamics.impl import SuperParticleSpawningDynamic
 
 
 @register_dynamic()
-class Seeding:
+class Seeding(SuperParticleSpawningDynamic):
     def __init__(
         self,
         *,
@@ -33,15 +34,10 @@ def register(self, builder):
         self.particulator = builder.particulator
 
     def post_register_setup_when_attributes_are_known(self):
-        if tuple(self.particulator.attributes.get_extensive_attribute_keys()) != tuple(
-            self.seeded_particle_extensive_attributes.keys()
-        ):
-            raise ValueError(
-                f"extensive attributes ({self.seeded_particle_extensive_attributes.keys()})"
-                " do not match those used in particulator"
-                f" ({self.particulator.attributes.get_extensive_attribute_keys()})"
-            )
-
+        SuperParticleSpawningDynamic.check_extensive_attribute_keys(
+            particulator_attributes=self.particulator.attributes,
+            spawned_attributes=self.seeded_particle_extensive_attributes,
+        )
         self.index = self.particulator.Index.identity_index(
             len(self.seeded_particle_multiplicity)
         )
@@ -86,9 +82,9 @@ def __call__(self):
                 # or if the number of super particles to inject
                 # is equal to the number of possible seeds
                 self.index.shuffle(self.u01)
-            self.particulator.seeding(
-                seeded_particle_index=self.index,
-                number_of_super_particles_to_inject=number_of_super_particles_to_inject,
-                seeded_particle_multiplicity=self.seeded_particle_multiplicity,
-                seeded_particle_extensive_attributes=self.seeded_particle_extensive_attributes,
+            self.particulator.spawn(
+                spawned_particle_index=self.index,
+                number_of_super_particles_to_spawn=number_of_super_particles_to_inject,
+                spawned_particle_multiplicity=self.seeded_particle_multiplicity,
+                spawned_particle_extensive_attributes=self.seeded_particle_extensive_attributes,
             )

PySDM/environments/__init__.py

@@ -8,3 +8,4 @@
 from .kinematic_1d import Kinematic1D
 from .kinematic_2d import Kinematic2D
 from .parcel import Parcel
+from .expansion_chamber import ExpansionChamber

PySDM/environments/expansion_chamber.py

@@ -0,0 +1,93 @@
+"""
+Zero-dimensional expansion chamber framework
+"""
+
+from typing import Optional, List
+
+from PySDM.environments.impl.moist_lagrangian import MoistLagrangian
+from PySDM.impl.mesh import Mesh
+from PySDM.environments.impl import register_environment
+
+
+@register_environment()
+class ExpansionChamber(MoistLagrangian):
+    def __init__(
+        self,
+        *,
+        dt,
+        volume: float,
+        initial_pressure: float,
+        initial_temperature: float,
+        initial_relative_humidity: float,
+        delta_pressure: float,
+        delta_time: float,
+        variables: Optional[List[str]] = None,
+        mixed_phase=False,
+    ):
+        variables = (variables or []) + ["rhod"]
+
+        super().__init__(dt, Mesh.mesh_0d(), variables, mixed_phase=mixed_phase)
+
+        self.dv = volume
+        self.initial_pressure = initial_pressure
+        self.initial_temperature = initial_temperature
+        self.initial_relative_humidity = initial_relative_humidity
+        self.delta_time = delta_time
+        self.dp_dt = delta_pressure / delta_time
+
+    def register(self, builder):
+        self.mesh.dv = self.dv
+
+        super().register(builder)
+
+        formulae = self.particulator.formulae
+        pv0 = (
+            self.initial_relative_humidity
+            * formulae.saturation_vapour_pressure.pvs_water(self.initial_temperature)
+        )
+        th_std = formulae.trivia.th_std(
+            p=self.initial_pressure, T=self.initial_temperature
+        )
+        initial_water_vapour_mixing_ratio = (
+            formulae.constants.eps * pv0 / (self.initial_pressure - pv0)
+        )
+
+        self["rhod"][:] = formulae.state_variable_triplet.rho_d(
+            p=self.initial_pressure,
+            water_vapour_mixing_ratio=initial_water_vapour_mixing_ratio,
+            theta_std=th_std,
+        )
+        self["thd"][:] = formulae.state_variable_triplet.th_dry(
+            th_std, initial_water_vapour_mixing_ratio
+        )
+        self["water_vapour_mixing_ratio"][:] = initial_water_vapour_mixing_ratio
+
+        self.post_register()
+
+    def advance_moist_vars(self):
+        """compute new values of dry density and thd, and write them to self._tmp and self["thd"]
+        assuming water-vapour mixing ratio is not altered by the expansion"""
+        dt = self.particulator.dt
+        t_new = (self.particulator.n_steps + 1) * dt
+        if t_new > self.delta_time:
+            return
+
+        formulae = self.particulator.formulae
+        p_new = self["p"][0] + self.dp_dt * dt
+        vapour_mixing_ratio = self["water_vapour_mixing_ratio"][0]
+        gg = (
+            1 - formulae.adiabatic_exponent.gamma(vapour_mixing_ratio)
+        ) / formulae.adiabatic_exponent.gamma(vapour_mixing_ratio)
+        T_new = self.initial_temperature * (self.initial_pressure / p_new) ** gg
+        wvmr_new = self._tmp["water_vapour_mixing_ratio"][
+            0
+        ]  # TODO #1492 - should _tmp or self[] be used?
+        pv_new = wvmr_new * p_new / (wvmr_new + formulae.constants.eps)
+        pd_new = p_new - pv_new
+
+        self._tmp["rhod"][:] = pd_new / T_new / formulae.constants.Rd
+        self["thd"][:] = formulae.trivia.th_std(p=pd_new, T=T_new)
+
+    def sync_moist_vars(self):
+        for var in self.variables:
+            self._tmp[var][:] = self[var][:]