homogeneous nucleation of liquid droplets & expansion chamber example

.gitignore

@@ -4,6 +4,9 @@ PySDM_examples/utils/temporary_files/*
 *.pkl
 *.vtu
 *.vts
+*.png
+*.csv
+*.xlsx
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

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
+        qv = self["water_vapour_mixing_ratio"][0]
+        gg = (
+            1 - formulae.adiabatic_exponent.gamma(qv)
+        ) / formulae.adiabatic_exponent.gamma(qv)
+        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][:]

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
 

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()

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.equilibrate_wet_radii 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,
         *,
@@ -60,7 +60,7 @@ def register(self, builder):
             rhod0, self.mass_of_dry_air
         )
 
-        Moist.register(self, builder)
+        super().register(builder)
 
         params = (
             self.initial_water_vapour_mixing_ratio,
@@ -74,11 +74,9 @@ def register(self, builder):
         self["rhod"][:] = params[2]
         self["z"][:] = params[3]
         self["t"][:] = params[4]
-
         self._tmp["water_vapour_mixing_ratio"][:] = params[0]
-        self.sync_parcel_vars()
-        Moist.sync(self)
-        self.notify()
+
+        self.post_register()
 
     def init_attributes(
         self,
@@ -108,7 +106,7 @@ def init_attributes(
             attributes["dry volume"] = dry_volume
         return attributes
 
-    def advance_parcel_vars(self):
+    def advance_moist_vars(self):
         dt = self.particulator.dt
         T = self["T"][0]
         p = self["p"][0]
@@ -143,21 +141,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()

PySDM/formulae.py

@@ -44,6 +44,7 @@ def __init__(  # pylint: disable=too-many-locals
         hydrostatics: str = "Default",
         freezing_temperature_spectrum: str = "Null",
         heterogeneous_ice_nucleation_rate: str = "Null",
+        homogeneous_liquid_nucleation_rate: str = "Null",
         fragmentation_function: str = "AlwaysN",
         isotope_equilibrium_fractionation_factors: str = "Null",
         isotope_meteoric_water_line: str = "Null",
@@ -57,6 +58,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
@@ -77,6 +79,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.fragmentation_function = fragmentation_function
         self.isotope_equilibrium_fractionation_factors = (
             isotope_equilibrium_fractionation_factors
@@ -90,6 +93,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

PySDM/physics/__init__.py

@@ -25,6 +25,7 @@
     fragmentation_function,
     freezing_temperature_spectrum,
     heterogeneous_ice_nucleation_rate,
+    homogeneous_liquid_nucleation_rate,
     hydrostatics,
     hygroscopicity,
     impl,
@@ -47,5 +48,6 @@
     air_dynamic_viscosity,
     terminal_velocity,
     bulk_phase_partitioning,
+    adiabatic_exponent,
 )
 from .constants import convert_to, in_unit, si

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

PySDM/physics/adiabatic_exponent/dry.py

@@ -0,0 +1,13 @@
+"""
+adiabatic exponent, dry air
+"""
+
+
+class Dry:
+    def __init__(self, _):
+        pass
+
+    # pylint: disable=too-many-arguments
+    @staticmethod
+    def gamma(const, qv):
+        return 1 + const.Rd / const.c_vd

PySDM/physics/adiabatic_exponent/moist_leading_terms.py

@@ -0,0 +1,18 @@
+"""
+adiabatic exponent, moist air, expanding to first order in qv, assuming qt=qv
+"""
+
+
+class MoistLeadingTerms:
+    def __init__(self, _):
+        pass
+
+    # pylint: disable=too-many-arguments
+    @staticmethod
+    def gamma(const, qv):
+        return (
+            1
+            + const.Rd / const.c_vd
+            + (const.Rv * qv) / (const.c_vd * (1 - qv))
+            - (const.Rd * qv * const.c_vv) / (const.c_vd**2 * (1 - qv))
+        )

PySDM/physics/constants_defaults.py

@@ -681,6 +681,9 @@ 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"]

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

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))

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