diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 75a6866658..0afe8144f2 100644
--- a/.github/workflows/tests.yml
+++ b/.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"
diff --git a/.gitignore b/.gitignore
index 4234ec2217..39fafa0be3 100644
--- a/.gitignore
+++ b/.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__/
diff --git a/PySDM/backends/impl_numba/methods/__init__.py b/PySDM/backends/impl_numba/methods/__init__.py
index 4c613462a2..b3b2573b9f 100644
--- a/PySDM/backends/impl_numba/methods/__init__.py
+++ b/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
diff --git a/PySDM/backends/impl_numba/methods/seeding_methods.py b/PySDM/backends/impl_numba/methods/spawning_methods.py
similarity index 54%
rename from PySDM/backends/impl_numba/methods/seeding_methods.py
rename to PySDM/backends/impl_numba/methods/spawning_methods.py
index 83fb64dc5a..653e27f191 100644
--- a/PySDM/backends/impl_numba/methods/seeding_methods.py
+++ b/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,
)
diff --git a/PySDM/backends/numba.py b/PySDM/backends/numba.py
index 930262a280..c522e23ebe 100644
--- a/PySDM/backends/numba.py
+++ b/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)
diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py
index d296bf4720..a8dc0b33cd 100644
--- a/PySDM/dynamics/__init__.py
+++ b/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
diff --git a/PySDM/dynamics/homogeneous_liquid_nucleation.py b/PySDM/dynamics/homogeneous_liquid_nucleation.py
new file mode 100644
index 0000000000..1f3869da61
--- /dev/null
+++ b/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
diff --git a/PySDM/dynamics/impl/__init__.py b/PySDM/dynamics/impl/__init__.py
index 2843142c69..455352bc1f 100644
--- a/PySDM/dynamics/impl/__init__.py
+++ b/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
diff --git a/PySDM/dynamics/impl/super_particle_spawning_dynamic.py b/PySDM/dynamics/impl/super_particle_spawning_dynamic.py
new file mode 100644
index 0000000000..64497d6541
--- /dev/null
+++ b/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()})"
+ )
diff --git a/PySDM/dynamics/seeding.py b/PySDM/dynamics/seeding.py
index 410d8abb20..894f783027 100644
--- a/PySDM/dynamics/seeding.py
+++ b/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,
)
diff --git a/PySDM/environments/__init__.py b/PySDM/environments/__init__.py
index f9e9d18f9c..e4fb48fdca 100644
--- a/PySDM/environments/__init__.py
+++ b/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
diff --git a/PySDM/environments/expansion_chamber.py b/PySDM/environments/expansion_chamber.py
new file mode 100644
index 0000000000..88ad00db27
--- /dev/null
+++ b/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][:]
diff --git a/PySDM/environments/impl/moist.py b/PySDM/environments/impl/moist.py
index bef7c02dc4..35009a166e 100644
--- a/PySDM/environments/impl/moist.py
+++ b/PySDM/environments/impl/moist.py
@@ -71,6 +71,8 @@ def _recalculate_temperature_pressure_relative_humidity(self, target):
)
def sync(self):
+ """fills the "predicted" set of thermodynamic variables with derived values freshly
+ computed from `self.get_thd()` and `self.get_water_vapour_mixing_ratio()`"""
target = self._tmp
target["water_vapour_mixing_ratio"].ravel(self.get_water_vapour_mixing_ratio())
target["thd"].ravel(self.get_thd())
@@ -108,6 +110,8 @@ def get_thd(self) -> np.ndarray:
raise NotImplementedError()
def notify(self):
+ """invalidates the "predicted" set of thermodynamic variables storing its
+ contents into the "current" set of vars"""
if self._values["predicted"] is None:
return
diff --git a/PySDM/environments/impl/moist_lagrangian.py b/PySDM/environments/impl/moist_lagrangian.py
new file mode 100644
index 0000000000..805d4b8d44
--- /dev/null
+++ b/PySDM/environments/impl/moist_lagrangian.py
@@ -0,0 +1,31 @@
+"""
+Zero-dimensional (Lagrangian) moist environment commons
+"""
+
+from PySDM.environments.impl.moist import Moist
+
+
+class MoistLagrangian(Moist):
+ """base class for moist Lagrangian environments (parcel, chamber, ...)"""
+
+ def get_thd(self):
+ return self["thd"]
+
+ def get_water_vapour_mixing_ratio(self):
+ return self["water_vapour_mixing_ratio"]
+
+ def sync(self):
+ self.sync_moist_vars()
+ self.advance_moist_vars()
+ super().sync()
+
+ def post_register(self):
+ self.sync_moist_vars()
+ super().sync()
+ self.notify()
+
+ def sync_moist_vars(self):
+ raise NotImplementedError()
+
+ def advance_moist_vars(self):
+ raise NotImplementedError()
diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py
index 1a501bcd45..9a31eb742e 100644
--- a/PySDM/environments/parcel.py
+++ b/PySDM/environments/parcel.py
@@ -6,7 +6,7 @@
import numpy as np
-from PySDM.environments.impl.moist import Moist
+from PySDM.environments.impl.moist_lagrangian import MoistLagrangian
from PySDM.impl.mesh import Mesh
from PySDM.initialisation.hygroscopic_equilibrium import (
default_rtol,
@@ -16,7 +16,7 @@
@register_environment()
-class Parcel(Moist): # pylint: disable=too-many-instance-attributes
+class Parcel(MoistLagrangian): # pylint: disable=too-many-instance-attributes
def __init__(
self,
*,
@@ -56,19 +56,16 @@ def register(self, builder):
rhod0, self.mass_of_dry_air
)
- Moist.register(self, builder)
+ super().register(builder)
self["water_vapour_mixing_ratio"][:] = self.initial_water_vapour_mixing_ratio
self["thd"][:] = formulae.trivia.th_std(pd0, self.T0)
self["rhod"][:] = rhod0
self["z"][:] = self.z0
-
self._tmp["water_vapour_mixing_ratio"][
:
] = self.initial_water_vapour_mixing_ratio
- self.sync_parcel_vars()
- Moist.sync(self)
- self.notify()
+ self.post_register()
def init_attributes(
self,
@@ -98,7 +95,7 @@ def init_attributes(
attributes["dry volume"] = dry_volume
return attributes
- def advance_parcel_vars(self):
+ def advance_moist_vars(self):
"""compute new values of displacement, dry-air density and volume,
and write them to self._tmp and self.mesh.dv"""
dt = self.particulator.dt
@@ -133,21 +130,10 @@ def advance_parcel_vars(self):
(self._tmp["rhod"][0] + self["rhod"][0]) / 2, self.mass_of_dry_air
)
- def get_thd(self):
- return self["thd"]
-
- def get_water_vapour_mixing_ratio(self):
- return self["water_vapour_mixing_ratio"]
-
- def sync_parcel_vars(self):
+ def sync_moist_vars(self):
self.delta_liquid_water_mixing_ratio = (
self._tmp["water_vapour_mixing_ratio"][0]
- self["water_vapour_mixing_ratio"][0]
)
for var in self.variables:
self._tmp[var][:] = self[var][:]
-
- def sync(self):
- self.sync_parcel_vars()
- self.advance_parcel_vars()
- super().sync()
diff --git a/PySDM/formulae.py b/PySDM/formulae.py
index 4b796f5784..df059a7655 100644
--- a/PySDM/formulae.py
+++ b/PySDM/formulae.py
@@ -47,6 +47,7 @@ def __init__( # pylint: disable=too-many-locals
hydrostatics: str = "ConstantGVapourMixingRatioAndThetaStd",
freezing_temperature_spectrum: str = "Null",
heterogeneous_ice_nucleation_rate: str = "Null",
+ homogeneous_liquid_nucleation_rate: str = "Null",
homogeneous_ice_nucleation_rate: str = "Null",
fragmentation_function: str = "AlwaysN",
isotope_equilibrium_fractionation_factors: str = "Null",
@@ -63,6 +64,7 @@ def __init__( # pylint: disable=too-many-locals
terminal_velocity: str = "GunnKinzer1949",
air_dynamic_viscosity: str = "ZografosEtAl1987",
bulk_phase_partitioning: str = "Null",
+ adiabatic_exponent: str = "Dry",
handle_all_breakups: bool = False,
):
# initialisation of the fields below is just to silence pylint and to enable code hints
@@ -87,6 +89,7 @@ def __init__( # pylint: disable=too-many-locals
self.hydrostatics = hydrostatics
self.freezing_temperature_spectrum = freezing_temperature_spectrum
self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate
+ self.homogeneous_liquid_nucleation_rate = homogeneous_liquid_nucleation_rate
self.homogeneous_ice_nucleation_rate = homogeneous_ice_nucleation_rate
self.fragmentation_function = fragmentation_function
self.isotope_equilibrium_fractionation_factors = (
@@ -105,6 +108,7 @@ def __init__( # pylint: disable=too-many-locals
self.air_dynamic_viscosity = air_dynamic_viscosity
self.terminal_velocity = terminal_velocity
self.bulk_phase_partitioning = bulk_phase_partitioning
+ self.adiabatic_exponent = adiabatic_exponent
self._components = tuple(
i
diff --git a/PySDM/particulator.py b/PySDM/particulator.py
index 5919feb805..dce9e2c6c5 100644
--- a/PySDM/particulator.py
+++ b/PySDM/particulator.py
@@ -445,39 +445,39 @@ def isotopic_fractionation(self, heavy_isotopes: tuple):
for isotope in heavy_isotopes:
self.attributes.mark_updated(f"moles_{isotope}")
- def seeding(
+ def spawn(
self,
*,
- seeded_particle_index,
- seeded_particle_multiplicity,
- seeded_particle_extensive_attributes,
- number_of_super_particles_to_inject,
+ spawned_particle_index,
+ spawned_particle_multiplicity,
+ spawned_particle_extensive_attributes,
+ number_of_super_particles_to_spawn,
):
n_null = self.n_sd - self.attributes.super_droplet_count
if n_null == 0:
raise ValueError(
- "No available seeds to inject. Please provide particles with nan filled attributes."
+ "No available null SDs to spawn. Please provide SDs with nan filled attributes."
)
- if number_of_super_particles_to_inject > n_null:
+ if number_of_super_particles_to_spawn > n_null:
raise ValueError(
- "Trying to inject more super particles than space available."
+ "Trying to spawn more super particles than space available."
)
- if number_of_super_particles_to_inject > len(seeded_particle_multiplicity):
+ if number_of_super_particles_to_spawn > len(spawned_particle_multiplicity):
raise ValueError(
- "Trying to inject multiple super particles with the same attributes. \
- Instead increase multiplicity of injected particles."
+ "Trying to spawn multiple super particles with the same attributes. \
+ Instead increase multiplicity of spawned particles."
)
- self.backend.seeding(
+ self.backend.spawning(
idx=self.attributes._ParticleAttributes__idx,
multiplicity=self.attributes["multiplicity"],
extensive_attributes=self.attributes.get_extensive_attribute_storage(),
- seeded_particle_index=seeded_particle_index,
- seeded_particle_multiplicity=seeded_particle_multiplicity,
- seeded_particle_extensive_attributes=seeded_particle_extensive_attributes,
- number_of_super_particles_to_inject=number_of_super_particles_to_inject,
+ spawned_particle_index=spawned_particle_index,
+ spawned_particle_multiplicity=spawned_particle_multiplicity,
+ spawned_particle_extensive_attributes=spawned_particle_extensive_attributes,
+ number_of_super_particles_to_spawn=number_of_super_particles_to_spawn,
)
self.attributes.reset_idx()
self.attributes.sanitize()
diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py
index f2e2bfb43f..22d6d11896 100644
--- a/PySDM/physics/__init__.py
+++ b/PySDM/physics/__init__.py
@@ -27,6 +27,7 @@
fragmentation_function,
freezing_temperature_spectrum,
heterogeneous_ice_nucleation_rate,
+ homogeneous_liquid_nucleation_rate,
homogeneous_ice_nucleation_rate,
hydrostatics,
hygroscopicity,
@@ -53,5 +54,6 @@
air_dynamic_viscosity,
terminal_velocity,
bulk_phase_partitioning,
+ adiabatic_exponent,
)
from .constants import convert_to, in_unit, si
diff --git a/PySDM/physics/adiabatic_exponent/__init__.py b/PySDM/physics/adiabatic_exponent/__init__.py
new file mode 100644
index 0000000000..d852b87f33
--- /dev/null
+++ b/PySDM/physics/adiabatic_exponent/__init__.py
@@ -0,0 +1,6 @@
+"""
+Alternative formulations of adiabatic exponent
+"""
+
+from .dry import Dry
+from .moist_leading_terms import MoistLeadingTerms
diff --git a/PySDM/physics/adiabatic_exponent/dry.py b/PySDM/physics/adiabatic_exponent/dry.py
new file mode 100644
index 0000000000..07db362c7f
--- /dev/null
+++ b/PySDM/physics/adiabatic_exponent/dry.py
@@ -0,0 +1,12 @@
+"""
+adiabatic exponent, dry air
+"""
+
+
+class Dry: # pylint: disable=too-few-public-methods
+ def __init__(self, _):
+ pass
+
+ @staticmethod
+ def gamma(const, qv): # pylint: disable=unused-argument
+ return 1 + const.Rd / const.c_vd
diff --git a/PySDM/physics/adiabatic_exponent/moist_leading_terms.py b/PySDM/physics/adiabatic_exponent/moist_leading_terms.py
new file mode 100644
index 0000000000..1e5e3bac23
--- /dev/null
+++ b/PySDM/physics/adiabatic_exponent/moist_leading_terms.py
@@ -0,0 +1,19 @@
+"""
+adiabatic exponent, moist air, expanding to first order in qv, assuming qt=qv
+"""
+
+
+class MoistLeadingTerms: # pylint: disable=too-few-public-methods
+ def __init__(self, _):
+ pass
+
+ @staticmethod
+ def gamma(const, vapour_mixing_ratio):
+ return (
+ 1
+ + const.Rd / const.c_vd
+ + (const.Rv * vapour_mixing_ratio)
+ / (const.c_vd * (1 - vapour_mixing_ratio))
+ - (const.Rd * vapour_mixing_ratio * const.c_vv)
+ / (const.c_vd**2 * (1 - vapour_mixing_ratio))
+ )
diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py
index 4144a04f5e..74f39c5eaa 100644
--- a/PySDM/physics/constants_defaults.py
+++ b/PySDM/physics/constants_defaults.py
@@ -359,6 +359,12 @@
J_HOM = np.nan
""" constant ice nucleation rates """
+J_LIQ_HOMO = np.nan
+""" constant liquid homogeneous-nucleation rate """
+
+R_LIQ_HOMO = np.nan
+""" constant particle size for liquid homogeneous nucleation """
+
STRAUB_E_D1 = 0.04 * si.cm
""" [Straub et al. 2010](https://doi.org/10.1175/2009JAS3175.1) """
STRAUB_MU2 = 0.095 * si.cm
@@ -770,8 +776,12 @@ def compute_derived_values(c: dict):
c["Rd_over_c_pd"] = c["Rd"] / c["c_pd"]
+ c["c_vd"] = c["c_pd"] - c["Rd"]
+ c["c_vv"] = c["c_pv"] - c["Rv"]
+
c["water_molar_volume"] = c["Mv"] / c["rho_w"]
c["rho_STP"] = c["p_STP"] / c["Rd"] / c["T_STP"]
c["H_u"] = c["M"] / c["p_STP"]
+ c["k_B"] = c["R_str"] / c["N_A"]
c["l_tri"] = c["L_tri"] / c["Mv"]
diff --git a/PySDM/physics/homogeneous_liquid_nucleation_rate/__init__.py b/PySDM/physics/homogeneous_liquid_nucleation_rate/__init__.py
new file mode 100644
index 0000000000..361803925f
--- /dev/null
+++ b/PySDM/physics/homogeneous_liquid_nucleation_rate/__init__.py
@@ -0,0 +1,7 @@
+"""
+homogeneous liquid droplet nucleation rate formulations
+"""
+
+from PySDM.impl.null_physics_class import Null
+from .cnt import CNT
+from .constant import Constant
diff --git a/PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py b/PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py
new file mode 100644
index 0000000000..23d47531f1
--- /dev/null
+++ b/PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py
@@ -0,0 +1,32 @@
+"""
+classical nucleation theory (CNT)
+formulae from [Seinfeld and Pandis](https://archive.org/details/0237-pdf-atmospheric-chemistry-and-physics-2nd-ed-j.-seinfeld-s.-pandis-wiley-2006-ww)
+Eq (11.47) and (11.52)
+""" # pylint: disable=line-too-long
+
+import numpy as np
+
+
+class CNT:
+ def __init__(self, _):
+ return
+
+ @staticmethod
+ def j_liq_homo(const, T, S, e_s):
+ mass_per_moleculue = const.Mv / const.N_A
+ volume_per_molecule = mass_per_moleculue / const.rho_w
+ N1 = (S * e_s) / (mass_per_moleculue * const.Rv * T) # molecules per m3
+ return (
+ ((2 * const.sgm_w) / (np.pi * mass_per_moleculue)) ** (1 / 2)
+ * (volume_per_molecule * N1**2 / S)
+ * np.exp(
+ (-16 * np.pi * volume_per_molecule**2 * const.sgm_w**3)
+ / (3 * const.k_B**3 * T**3 * np.log(S) ** 2)
+ )
+ )
+
+ @staticmethod
+ def r_liq_homo(const, T, S):
+ mass_per_moleculue = const.Mv / const.N_A
+ volume_per_molecule = mass_per_moleculue / const.rho_w
+ return (2 * const.sgm_w * volume_per_molecule) / (const.k_B * T * np.log(S))
diff --git a/PySDM/physics/homogeneous_liquid_nucleation_rate/constant.py b/PySDM/physics/homogeneous_liquid_nucleation_rate/constant.py
new file mode 100644
index 0000000000..49dd74796f
--- /dev/null
+++ b/PySDM/physics/homogeneous_liquid_nucleation_rate/constant.py
@@ -0,0 +1,19 @@
+"""
+constant rate formulation (for tests)
+"""
+
+import numpy as np
+
+
+class Constant:
+ def __init__(self, const):
+ assert np.isfinite(const.J_LIQ_HOMO)
+ assert np.isfinite(const.R_LIQ_HOMO)
+
+ @staticmethod
+ def j_liq_homo(const, T, S, e_s): # pylint: disable=unused-argument
+ return const.J_LIQ_HOMO
+
+ @staticmethod
+ def r_liq_homo(const, T, S): # pylint: disable=unused-argument
+ return const.R_LIQ_HOMO
diff --git a/PySDM/physics/impl/fake_unit_registry.py b/PySDM/physics/impl/fake_unit_registry.py
index ddac31daef..1665cd0bb9 100644
--- a/PySDM/physics/impl/fake_unit_registry.py
+++ b/PySDM/physics/impl/fake_unit_registry.py
@@ -30,6 +30,7 @@ def __init__(self, si):
"hour",
"newton",
"watt",
+ "dyne",
):
self.__setattr__(prefix + unit, _fake(si.__getattr__(prefix + unit)))
self.__setattr__(
@@ -54,5 +55,6 @@ def __init__(self, si):
"bar", # note: "b" is barn !!!
"N",
"W",
+ "dyn",
):
self.__setattr__(prefix + unit, _fake(si.__getattr__(prefix + unit)))
diff --git a/docs/bibliography.json b/docs/bibliography.json
index 60eb54447e..a54f958574 100644
--- a/docs/bibliography.json
+++ b/docs/bibliography.json
@@ -730,6 +730,23 @@
"title": "Vapor Pressure Formulation for Water in Range 0 to 100 °C. A Revision",
"label": "Wexler 1976 (J. Res. NBS A Phys. Ch. 80A)"
},
+ "https://doi.org/10.48550/arXiv.2501.01467": {
+ "usages": [
+ "examples/PySDM_examples/Erinin_et_al_2025/__init__.py",
+ "examples/PySDM_examples/Erinin_et_al_2025/aerosol.py",
+ "tests/smoke_tests/chamber/erinin_et_al_2025/test_fig_3.py",
+ "examples/PySDM_examples/Erinin_et_al_2025/expansion_experiment.ipynb"],
+ "title": "Droplet Nucleation In a Rapid Expansion Aerosol Chamber",
+ "label": "Erinin et al. 2024 (arXiv:2501.01467)"
+ },
+ "https://archive.org/details/0237-pdf-atmospheric-chemistry-and-physics-2nd-ed-j.-seinfeld-s.-pandis-wiley-2006-ww": {
+ "usages": [
+ "PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py",
+ "PySDM/physics/constants_defaults.py"
+ ],
+ "title": "Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd Edition",
+ "label": "Seinfeld & Pandis 2006 (Wiley)"
+ },
"https://doi.org/10.1080/10789669.2008.10391032": {
"usages": ["PySDM/physics/constants_defaults.py"],
"title": "A Twenty-First Century Molar Mass for Dry Air",
@@ -784,11 +801,6 @@
"title": "Microphysics of Clouds and Precipitation",
"label": "Pruppacher & Klett 2010 (Springer Dordrecht)"
},
- "https://archive.org/details/0237-pdf-atmospheric-chemistry-and-physics-2nd-ed-j.-seinfeld-s.-pandis-wiley-2006-ww": {
- "usages": ["PySDM/physics/constants_defaults.py"],
- "title": "Microphysics of Clouds and Precipitation",
- "label": "Seinfeld & Pandis 2006 (Wiley)"
- },
"https://doi.org/10.1038/187857a0": {
"usages": [
"PySDM/physics/constants_defaults.py",
diff --git a/examples/PySDM_examples/Erinin_et_al_2025/__init__.py b/examples/PySDM_examples/Erinin_et_al_2025/__init__.py
new file mode 100644
index 0000000000..3916bea269
--- /dev/null
+++ b/examples/PySDM_examples/Erinin_et_al_2025/__init__.py
@@ -0,0 +1,4 @@
+# pylint: disable=invalid-name
+"""
+expansion chamber example based on [Erinin et al. 2025](https://doi.org/10.48550/arXiv.2501.01467)
+"""
diff --git a/examples/PySDM_examples/Erinin_et_al_2025/aerosol.py b/examples/PySDM_examples/Erinin_et_al_2025/aerosol.py
new file mode 100644
index 0000000000..36dadaa51d
--- /dev/null
+++ b/examples/PySDM_examples/Erinin_et_al_2025/aerosol.py
@@ -0,0 +1,47 @@
+"""aerosol parameters from [Erinin et al. 2025](https://doi.org/10.48550/arXiv.2501.01467)"""
+
+from chempy import Substance
+from pystrict import strict
+
+from PySDM.initialisation import spectra
+from PySDM.initialisation.aerosol_composition import DryAerosolMixture
+from PySDM.physics import si
+
+
+@strict
+class AerosolChamber(DryAerosolMixture):
+ def __init__(
+ self,
+ water_molar_volume: float,
+ N: float = 1e5 / si.cm**3,
+ ):
+ mode1 = {"CaCO3": 1.0}
+
+ super().__init__(
+ compounds=("CaCO3",),
+ molar_masses={
+ "CaCO3": Substance.from_formula("CaCO3").mass * si.g / si.mole
+ },
+ densities={
+ "CaCO3": 2.71 * si.g / si.cm**3,
+ },
+ is_soluble={
+ "CaCO3": True,
+ },
+ ionic_dissociation_phi={
+ "CaCO3": 1,
+ },
+ )
+ self.modes = (
+ {
+ "kappa": self.kappa(
+ mass_fractions=mode1,
+ water_molar_volume=water_molar_volume,
+ ),
+ "spectrum": spectra.Gaussian(
+ norm_factor=N,
+ loc=316 / 2 * si.nm,
+ scale=257 / 2 * si.nm,
+ ),
+ },
+ )
diff --git a/examples/PySDM_examples/Erinin_et_al_2025/expansion_experiment.ipynb b/examples/PySDM_examples/Erinin_et_al_2025/expansion_experiment.ipynb
new file mode 100644
index 0000000000..d4b936771a
--- /dev/null
+++ b/examples/PySDM_examples/Erinin_et_al_2025/expansion_experiment.ipynb
@@ -0,0 +1,10586 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "[](https://github.com/open-atmos/PySDM/blob/main/examples/PySDM_examples/Erinin_et_al_2025/expansion_experiment.ipynb)\n",
+ "[](https://mybinder.org/v2/gh/open-atmos/PySDM.git/main?urlpath=lab/tree/examples/PySDM_examples/Erinin_et_al_2025/expansion_experiment.ipynb)\n",
+ "[](https://colab.research.google.com/github/open-atmos/PySDM/blob/main/examples/PySDM_examples/Erinin_et_al_2025/expansion_experiment.ipynb)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "modeling expansion-chamber experiments described in [Erinin et al. 2025](https://doi.org/10.48550/arXiv.2501.01467)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-30T20:21:31.584604Z",
+ "start_time": "2025-01-30T20:21:31.573353Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "if 'google.colab' in sys.modules:\n",
+ " !pip --quiet install open-atmos-jupyter-utils\n",
+ " from open_atmos_jupyter_utils import pip_install_on_colab\n",
+ " pip_install_on_colab('PySDM-examples')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-30T20:21:37.280272Z",
+ "start_time": "2025-01-30T20:21:31.850352Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from matplotlib import pyplot\n",
+ "from open_atmos_jupyter_utils import show_plot\n",
+ "from PySDM import Formulae\n",
+ "from PySDM.physics import si, in_unit\n",
+ "from PySDM import products as PySDM_products\n",
+ "from PySDM_examples.Erinin_et_al_2025.aerosol import AerosolChamber\n",
+ "from PySDM_examples.Erinin_et_al_2025.expansion_simulation import run_expansion"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dry_radius_bin_edges = np.geomspace(50 * si.nm, 2000 * si.nm, 40, endpoint=False)\n",
+ "wet_radius_bin_edges = np.geomspace(1 * si.um, 40 * si.um, 40, endpoint=False)\n",
+ "\n",
+ "products = (\n",
+ " PySDM_products.SuperDropletCountPerGridbox(name=\"SD count\"),\n",
+ " PySDM_products.WaterMixingRatio(unit=\"g/kg\", name=r\"$q_\\ell$\"),\n",
+ " PySDM_products.PeakSupersaturation(name=\"Saturation ratio\"),\n",
+ " PySDM_products.AmbientRelativeHumidity(name=\"RH\"),\n",
+ " PySDM_products.AmbientTemperature(name=\"Temperature\", var=\"T\"),\n",
+ " PySDM_products.AmbientPressure(name=\"Pressure\", var=\"p\", unit=\"hPa\"),\n",
+ " PySDM_products.AmbientWaterVapourMixingRatio(\n",
+ " unit=\"g/kg\", name=\"$q_v$\", var=\"water_vapour_mixing_ratio\"\n",
+ " ),\n",
+ " PySDM_products.Time(name=\"t\"),\n",
+ " PySDM_products.ParticleSizeSpectrumPerVolume(\n",
+ " name=\"dry:dN/dR\",\n",
+ " unit=\"m^-3 m^-1\",\n",
+ " radius_bins_edges=dry_radius_bin_edges,\n",
+ " dry=True,\n",
+ " ),\n",
+ " PySDM_products.ParticleSizeSpectrumPerVolume(\n",
+ " name=\"wet:dN/dR\",\n",
+ " unit=\"m^-3 m^-1\",\n",
+ " radius_bins_edges=wet_radius_bin_edges,\n",
+ " dry=False,\n",
+ " ),\n",
+ " PySDM_products.ActivatedEffectiveRadius(\n",
+ " name=\"act_reff\", unit=\"um\", count_activated=True, count_unactivated=False\n",
+ " ),\n",
+ " PySDM_products.EffectiveRadius(\n",
+ " name=\"$r_{eff}$\",\n",
+ " unit=\"um\",\n",
+ " ),\n",
+ " PySDM_products.ParticleConcentration(\n",
+ " name=\"$N_d$\",\n",
+ " unit=\"cm^-3\",\n",
+ " # radius_range=(0.5 * si.um, 25 * si.um),\n",
+ " ),\n",
+ ") "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "FORMULAE = Formulae(\n",
+ " adiabatic_exponent=\"MoistLeadingTerms\",\n",
+ " homogeneous_liquid_nucleation_rate='CNT',\n",
+ ")\n",
+ "CONST = FORMULAE.constants\n",
+ "\n",
+ "N_SD_PER_MODE = 20\n",
+ "DT = .05 * si.s\n",
+ "\n",
+ "C0_arr = np.array([1e2, 2.5e2, 5e2, 7.5e2, 1e3]) / si.cm**3\n",
+ "C0_arr = np.array([4e2, 4.25e2, 4.5e2, 5e2, 5.5e2, 6e2]) / si.cm**3"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-30T20:23:04.715504Z",
+ "start_time": "2025-01-30T20:21:37.309294Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "C0=4.00e+08\n",
+ "C0=4.25e+08\n",
+ "C0=4.50e+08\n",
+ "C0=5.00e+08\n",
+ "C0=5.50e+08\n",
+ "C0=6.00e+08\n"
+ ]
+ }
+ ],
+ "source": [
+ "output_all = {}\n",
+ "\n",
+ "for Na in C0_arr:\n",
+ " key = f\"C0={Na:.2e}\"\n",
+ " print(key)\n",
+ " aerosol = AerosolChamber(\n",
+ " water_molar_volume=CONST.Mv / CONST.rho_w,\n",
+ " N=Na,\n",
+ " )\n",
+ " output = run_expansion(\n",
+ " formulae=FORMULAE,\n",
+ " aerosol=aerosol,\n",
+ " n_sd_per_mode=N_SD_PER_MODE,\n",
+ " n_sd_homo_liq_nucleation=100,\n",
+ " total_time=4*si.s,\n",
+ " dt=DT,\n",
+ " products=products,\n",
+ " )\n",
+ "\n",
+ " output_all[key] = output\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "4.021660843931165\n",
+ "0.33964378383596744\n",
+ "0.0197293310945584\n",
+ "7.870157403148062e-05\n",
+ "4.28580000171432e-07\n",
+ "0.0\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "d5c1bb95a07e4ff7b732ae5711345dbf",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./C0_sweep_traces.pdf
\"), HT…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "variables = [\"Pressure\", \"Temperature\", \"Saturation ratio\", \"SD count\", r\"$q_\\ell$\", \"$q_v$\", \"$r_{eff}$\", \"$N_d$\"]\n",
+ "fig, axes = pyplot.subplots(2,4,figsize=(12,5), sharex=True, sharey=False, constrained_layout=True)\n",
+ "\n",
+ "for Na in C0_arr:\n",
+ " key = f\"C0={Na:.2e}\"\n",
+ " output = output_all[key]\n",
+ " init_mult = output.attributes[\"multiplicity\"][0]\n",
+ " final_mult = output.attributes[\"multiplicity\"][-1]\n",
+ " print((np.sum(final_mult) - np.sum(init_mult)) / np.sum(init_mult))\n",
+ " for i, ax in enumerate(axes.flatten()):\n",
+ " y = np.array(output.profile[variables[i]])\n",
+ " ax.plot(output.profile[\"t\"], y, label=f\"{in_unit(Na, 1/si.cm**3):.2e}\")\n",
+ " if i == 7:\n",
+ " ax.set_yscale(\"log\")\n",
+ " ax.set_ylim(3e2,8e2)\n",
+ " ax.set_xlabel(f\"Time [{output.units['t']}]\")\n",
+ " ax.set_ylabel(f\"{variables[i]} [{output.units[variables[i]]}]\")\n",
+ "\n",
+ "axes[0, 0].legend(title=\"$C_0$ [cm$^{-3}$]\")\n",
+ "show_plot(\"C0_sweep_traces.pdf\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-30T20:23:04.785100072Z",
+ "start_time": "2025-01-26T23:04:07.407334Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "dpi=100.0 case='base'\n",
+ "dpi=100.0 case='clean'\n",
+ "dpi=100.0 case='dry'\n",
+ "dpi=10000.0 case='base'\n",
+ "dpi=10000.0 case='clean'\n",
+ "dpi=10000.0 case='dry'\n",
+ "dpi=20000.0 case='base'\n",
+ "dpi=20000.0 case='clean'\n",
+ "dpi=20000.0 case='dry'\n",
+ "dpi=30000.0 case='base'\n",
+ "dpi=30000.0 case='clean'\n",
+ "dpi=30000.0 case='dry'\n",
+ "dpi=40000.0 case='base'\n",
+ "dpi=40000.0 case='clean'\n",
+ "dpi=40000.0 case='dry'\n",
+ "dpi=50000.0 case='base'\n",
+ "dpi=50000.0 case='clean'\n",
+ "dpi=50000.0 case='dry'\n"
+ ]
+ }
+ ],
+ "source": [
+ "DP = np.insert(np.linspace(100, 500, 5), 0, 1) * si.hPa\n",
+ "P0 = 1000 * si.hPa\n",
+ "\n",
+ "cases = {\n",
+ " 'base': {'N': 2000 / si.cm**3, 'RH0': 0.5, 'color': 'green'},\n",
+ " 'clean': {'N': 1 / si.cm**3, 'RH0': 0.5, 'color': 'turquoise'},\n",
+ " 'dry': {'N': 1 / si.cm**3, 'RH0': 0, 'color': 'black'},\n",
+ "}\n",
+ "output = {\n",
+ " 'Smax': {case: [] for case in cases},\n",
+ " 'Tmin': {case: [] for case in cases},\n",
+ "}\n",
+ "for dpi in DP:\n",
+ " for case, case_params in cases.items():\n",
+ " print(f\"{dpi=} {case=}\")\n",
+ " aerosol = AerosolChamber(water_molar_volume=CONST.Mv / CONST.rho_w, N=case_params['N'])\n",
+ " out = run_expansion(\n",
+ " formulae=FORMULAE,\n",
+ " aerosol=aerosol,\n",
+ " n_sd_per_mode=N_SD_PER_MODE,\n",
+ " RH0=case_params['RH0'],\n",
+ " pf=(P0 - dpi),\n",
+ " dt=DT,\n",
+ " p0=P0,\n",
+ " products=products\n",
+ " )\n",
+ " output['Tmin'][case].append(np.nanmin(out.profile[\"Temperature\"]))\n",
+ " output['Smax'][case].append(np.nanmax(out.profile[\"Saturation ratio\"]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-30T20:18:02.702242840Z",
+ "start_time": "2025-01-26T23:10:40.873425Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "def plot(_axes, _output, _cases, var):\n",
+ " for _case, _case_params in _cases.items():\n",
+ " _axes.plot(\n",
+ " DP / P0,\n",
+ " _output[var][_case],\n",
+ " color=_case_params['color'],\n",
+ " marker=\"o\",\n",
+ " ls=\"--\",\n",
+ " label=f\"$S_0={_case_params['RH0']}, C_0={in_unit(_case_params['N'], si.cm**-3)}$ cm$^{{-3}}$\"\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 35,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "ee6aed70075c4f35ae142be20468587b",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./Tmin_dp.pdf
\"), HTML(value=\"\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "4970459fb58448c9bdd12c1471dea3d0",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./dp_sweep.pdf
\"), HTML(value=\"\n\n\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "b9ff2d6977b84d08b84d5c3dabb2139e",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HTML(value=\"./PF_sweep_traces.pdf
\")"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "variables = [\"Pressure\", \"Temperature\", \"Saturation ratio\", \"SD count\", r\"$q_\\ell$\", \"$q_v$\", \"$r_{eff}$\", \"$N_d$\"]\n",
+ "fig, axes = pyplot.subplots(2,4,figsize=(12,5), sharex=True, sharey=False, constrained_layout=True)\n",
+ "\n",
+ "for PFi in PF_arr:\n",
+ " key = f\"pf={PFi:.2e}\"\n",
+ " output = output_all[key]\n",
+ " init_mult = output.attributes[\"multiplicity\"][0]\n",
+ " final_mult = output.attributes[\"multiplicity\"][-1]\n",
+ " print((np.sum(final_mult) - np.sum(init_mult)) / np.sum(init_mult))\n",
+ " for i, ax in enumerate(axes.flatten()):\n",
+ " y = np.array(output.profile[variables[i]])\n",
+ " ax.plot(output.profile[\"t\"], y, label=f\"{in_unit(PFi, si.hPa):.2e}\")\n",
+ " if i == 7:\n",
+ " ax.set_yscale(\"log\")\n",
+ " # ax.set_ylim(3e2,8e2)\n",
+ " ax.set_xlabel(f\"Time [{output.units['t']}]\")\n",
+ " ax.set_ylabel(f\"{variables[i]} [{output.units[variables[i]]}]\")\n",
+ "\n",
+ "axes[0, 0].legend(title=\"$p_f$ [hPa]\")\n",
+ "fig.set_facecolor('white')\n",
+ "show_plot(\"PF_sweep_traces.pdf\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.9"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "b14f34a08619f4a218d80d7380beed3f0c712c89ff93e7183219752d640ed427"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/PySDM_examples/Erinin_et_al_2025/expansion_simulation.py b/examples/PySDM_examples/Erinin_et_al_2025/expansion_simulation.py
new file mode 100644
index 0000000000..fd4409c8a5
--- /dev/null
+++ b/examples/PySDM_examples/Erinin_et_al_2025/expansion_simulation.py
@@ -0,0 +1,151 @@
+from collections import namedtuple
+
+import numpy as np
+from PySDM import Builder
+from PySDM.backends import CPU
+from PySDM.dynamics import (
+ AmbientThermodynamics,
+ Condensation,
+ HomogeneousLiquidNucleation,
+)
+from PySDM.environments import ExpansionChamber
+from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii
+from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
+from PySDM.physics import si
+
+
+def run_expansion(
+ *,
+ formulae,
+ aerosol,
+ n_sd_per_mode,
+ n_sd_homo_liq_nucleation=100,
+ RH0=0.7,
+ T0=296 * si.K,
+ p0=1000 * si.hPa,
+ pf=500 * si.hPa,
+ delta_time=2 * si.s,
+ total_time=4 * si.s,
+ dt=0.1 * si.s,
+ volume=0.14 * si.m**3,
+ products=None,
+):
+ n_steps = int(np.ceil(total_time / dt))
+
+ env = ExpansionChamber(
+ dt=dt,
+ initial_pressure=p0,
+ delta_pressure=pf - p0,
+ initial_temperature=T0,
+ initial_relative_humidity=RH0,
+ delta_time=delta_time,
+ volume=volume,
+ )
+
+ n_sd = n_sd_per_mode * len(aerosol.modes) + n_sd_homo_liq_nucleation
+
+ builder = Builder(
+ backend=CPU(formulae, override_jit_flags={"parallel": False}),
+ n_sd=n_sd,
+ environment=env,
+ )
+ builder.add_dynamic(AmbientThermodynamics())
+ builder.add_dynamic(Condensation(adaptive=True))
+ builder.add_dynamic(HomogeneousLiquidNucleation())
+ builder.request_attribute("critical saturation")
+
+ attributes = {
+ k: np.empty(0) for k in ("dry volume", "kappa times dry volume", "multiplicity")
+ }
+ for mode in aerosol.modes:
+ kappa, spectrum = mode["kappa"]["Constant"], mode["spectrum"]
+ r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd_per_mode)
+ v_dry = builder.formulae.trivia.volume(radius=r_dry)
+ attributes["multiplicity"] = np.append(
+ attributes["multiplicity"],
+ concentration * builder.particulator.environment.dv,
+ )
+ attributes["dry volume"] = np.append(attributes["dry volume"], v_dry)
+ attributes["kappa times dry volume"] = np.append(
+ attributes["kappa times dry volume"], v_dry * kappa
+ )
+
+ r_wet = equilibrate_wet_radii(
+ r_dry=builder.formulae.trivia.radius(volume=attributes["dry volume"]),
+ environment=builder.particulator.environment,
+ kappa_times_dry_volume=attributes["kappa times dry volume"],
+ )
+ attributes["volume"] = builder.formulae.trivia.volume(radius=r_wet)
+
+ particulator = builder.build(
+ attributes={
+ k: np.pad(
+ array=v,
+ pad_width=(0, n_sd_homo_liq_nucleation),
+ mode="constant",
+ constant_values=np.nan if k == "multiplicity" else 0,
+ )
+ for k, v in attributes.items()
+ },
+ products=products or (),
+ )
+
+ output = {product.name: [] for product in particulator.products.values()}
+ output_attributes = {
+ k: []
+ for k in (
+ "multiplicity",
+ "volume",
+ "critical volume",
+ "critical saturation",
+ )
+ }
+
+ for _ in range(n_steps):
+ particulator.run(steps=1)
+ for product in particulator.products.values():
+ if product.name in ("dry:dN/dR", "wet:dN/dR"):
+ continue
+ value = product.get()
+ if product.name == "t":
+ output[product.name].append(value)
+ else:
+ output[product.name].append(value[0])
+ mult = particulator.attributes["multiplicity"].to_ndarray(raw=True)
+ for key, attr in output_attributes.items():
+ if key == "multiplicity":
+ attr.append(mult)
+ continue
+ data = particulator.attributes[key].to_ndarray(raw=True)
+ data[mult == 0] = np.nan
+ attr.append(data)
+
+ dry_spectrum = particulator.products["dry:dN/dR"].get()
+ wet_spectrum = particulator.products["wet:dN/dR"].get()
+
+ Output = namedtuple(
+ "Output",
+ (
+ "profile",
+ "attributes",
+ "aerosol",
+ "dry_spectrum",
+ "wet_spectrum",
+ "units",
+ ),
+ )
+ return Output(
+ profile=output,
+ attributes=output_attributes,
+ aerosol=aerosol,
+ dry_spectrum=dry_spectrum,
+ wet_spectrum=wet_spectrum,
+ units={
+ key: (
+ product.unit[2:]
+ if product.unit[0:2] == "1 "
+ else product.unit[4:] if product.unit[0:4] == "1.0 " else product.unit
+ ).replace("**", "^")
+ for key, product in particulator.products.items()
+ },
+ )
diff --git a/examples/PySDM_examples/Erinin_et_al_2025/fig_3.ipynb b/examples/PySDM_examples/Erinin_et_al_2025/fig_3.ipynb
new file mode 100644
index 0000000000..49e5526837
--- /dev/null
+++ b/examples/PySDM_examples/Erinin_et_al_2025/fig_3.ipynb
@@ -0,0 +1,1745 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "f523b05e-2bd0-4454-8256-6d7d4e9fbfc8",
+ "metadata": {},
+ "source": [
+ "[](https://github.com/open-atmos/PySDM/blob/main/examples/PySDM_examples/Erinin_et_al_2025/fig_3.ipynb)\n",
+ "[](https://mybinder.org/v2/gh/open-atmos/PySDM.git/main?urlpath=lab/tree/examples/PySDM_examples/Erinin_et_al_2025/fig_3.ipynb)\n",
+ "[](https://colab.research.google.com/github/open-atmos/PySDM/blob/main/examples/PySDM_examples/Erinin_et_al_2025/fig_3.ipynb)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "13f49050-5018-4897-af86-af897858d521",
+ "metadata": {},
+ "source": [
+ "TODO #1492"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "c84406ae-5606-4ea1-888c-be41826f6a9d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "if 'google.colab' in sys.modules:\n",
+ " !pip --quiet install open-atmos-jupyter-utils\n",
+ " from open_atmos_jupyter_utils import pip_install_on_colab\n",
+ " pip_install_on_colab('PySDM-examples')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "initial_id",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-25T23:31:05.929440Z",
+ "start_time": "2025-01-25T23:31:02.733943Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from matplotlib import pyplot\n",
+ "import numpy as np\n",
+ "from open_atmos_jupyter_utils import show_plot\n",
+ "from PySDM.environments import ExpansionChamber\n",
+ "from PySDM import products\n",
+ "from PySDM import Builder, Formulae\n",
+ "from PySDM.physics.constants import si, in_unit, PER_CENT\n",
+ "from PySDM.dynamics import AmbientThermodynamics\n",
+ "from PySDM.backends import CPU\n",
+ "from PySDM_examples.utils import BasicSimulation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "7e6e3c6c48db80ee",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-26T23:05:36.133638Z",
+ "start_time": "2025-01-26T23:05:36.125604Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "class Simulation(BasicSimulation): # pylint: disable=too-few-public-methods\n",
+ " \"\"\" minimal simulation setup \"\"\"\n",
+ " def __init__(self):\n",
+ " builder = Builder(n_sd=0, backend=CPU(), environment=ExpansionChamber(\n",
+ " dt = .1 * si.s,\n",
+ " volume=.14 * si.m**3,\n",
+ " initial_pressure=1 * si.bar,\n",
+ " initial_temperature=293.15 * si.K,\n",
+ " initial_relative_humidity=0.5,\n",
+ " delta_time=1 * si.s,\n",
+ " delta_pressure=-.54 * si.bar,\n",
+ " ))\n",
+ " builder.add_dynamic(AmbientThermodynamics())\n",
+ " super().__init__(builder.build(\n",
+ " attributes={\n",
+ " \"multiplicity\": np.array([1]),\n",
+ " \"water mass\": np.array([1 * si.ug]),\n",
+ " },\n",
+ " products=(\n",
+ " products.AmbientPressure(var='p'),\n",
+ " products.AmbientTemperature(var='T'),\n",
+ " products.AmbientRelativeHumidity(var='RH'),\n",
+ " products.Time(),\n",
+ " )\n",
+ " ))\n",
+ "\n",
+ " def run(self):\n",
+ " return self._run(nt=int(5 * si.s / self.particulator.environment.dt), steps_per_output_interval=1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "d658455c740ae2de",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-26T23:05:45.266151Z",
+ "start_time": "2025-01-26T23:05:37.527519Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "simulation = Simulation()\n",
+ "output = simulation.run()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "d94f32bcd47f8205",
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-01-26T23:09:45.305277Z",
+ "start_time": "2025-01-26T23:09:42.691896Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "f7bc96feb07a4b7a9a6738491c8d97c6",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./fig_3.pdf
\"), HTML(value=\"