new example: Toon et al. 1980 (hello-world for using PySDM for extraterrestrial clouds, here for Titan); variable-g/constant-T hydrostatics formulae (#142)

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

Author: 3 people

.github/workflows/tests.yml

@@ -213,7 +213,7 @@ jobs:
       matrix:
         platform: [ubuntu-24.04, macos-13, windows-latest]
         python-version: ["3.9", "3.12"]
-        test-suite: [ "chemistry_freezing", "isotopes", "condensation_a", "condensation_b", "condensation_c", "coagulation", "breakup", "multi-process_a", "multi-process_b", "multi-process_c", "multi-process_d", "multi-process_e"]
+        test-suite: [ "chemistry_freezing_extraterrestrial", "isotopes", "condensation_a", "condensation_b", "condensation_c", "coagulation", "breakup", "multi-process_a", "multi-process_b", "multi-process_c", "multi-process_d", "multi-process_e"]
       fail-fast: false
     runs-on: ${{ matrix.platform }}
     timeout-minutes: 35

PySDM/environments/parcel.py

@@ -114,7 +114,6 @@ def advance_parcel_vars(self):
 
         # derivative evaluated at p_old, T_old, mixrat_mid, w_mid
         drho_dz = formulae.hydrostatics.drho_dz(
-            g=formulae.constants.g_std,
             p=p,
             T=T,
             water_vapour_mixing_ratio=water_vapour_mixing_ratio,

PySDM/formulae.py

@@ -41,7 +41,7 @@ def __init__(  # pylint: disable=too-many-locals
         ventilation: str = "Neglect",
         state_variable_triplet: str = "LibcloudphPlusPlus",
         particle_advection: str = "ImplicitInSpace",
-        hydrostatics: str = "Default",
+        hydrostatics: str = "ConstantGVapourMixingRatioAndThetaStd",
         freezing_temperature_spectrum: str = "Null",
         heterogeneous_ice_nucleation_rate: str = "Null",
         fragmentation_function: str = "AlwaysN",

PySDM/physics/constants_defaults.py

@@ -606,6 +606,9 @@
 g_std = sci.g * si.metre / si.second**2
 """ standard gravitational acceleration (value from SciPy) """
 
+celestial_body_radius = np.nan
+""" radius of the considered celestial body (Earth, Titan, ...) """
+
 bulk_phase_partitioning_T_cold = 235 * si.K
 """ [Kaul et al. 2015](https://doi.org/10.1175/MWR-D-14-00319.1) """
 bulk_phase_partitioning_T_warm = 273 * si.K

PySDM/physics/hydrostatics/__init__.py

@@ -2,4 +2,7 @@
 Helper routines for hydrostatic pressure and density profiles
 """
 
-from .default import Default
+from .constant_g_vapour_mixing_ratio_and_theta_std import (
+    ConstantGVapourMixingRatioAndThetaStd,
+)
+from .variable_g_isothermal import VariableGIsothermal

PySDM/physics/hydrostatics/default.py …t_g_vapour_mixing_ratio_and_theta_std.pyPySDM/physics/hydrostatics/default.py renamed to PySDM/physics/hydrostatics/constant_g_vapour_mixing_ratio_and_theta_std.py PySDM/physics/hydrostatics/default.py renamed to PySDM/physics/hydrostatics/constant_g_vapour_mixing_ratio_and_theta_std.py

@@ -1,18 +1,19 @@
 """
-# TODO #407
+hydrostatic profile assuming constant g, constant water vapour mixing ratio
+and constant (standard) potential temperature
 """
 
 import numpy as np
 
 
-class Default:
+class ConstantGVapourMixingRatioAndThetaStd:
     def __init__(self, _):
         pass
 
     # pylint: disable=too-many-arguments
     @staticmethod
     def drho_dz(
-        const, g, p, T, water_vapour_mixing_ratio, lv, d_liquid_water_mixing_ratio__dz=0
+        const, p, T, water_vapour_mixing_ratio, lv, d_liquid_water_mixing_ratio__dz=0
     ):
         Rq = const.Rv / (1 / water_vapour_mixing_ratio + 1) + const.Rd / (
             1 + water_vapour_mixing_ratio
@@ -22,21 +23,21 @@ def drho_dz(
         )
         rho = p / Rq / T
         return (
-            g / T * rho * (Rq / cp - 1)
+            const.g_std / T * rho * (Rq / cp - 1)
             - p * lv / cp / T**2 * d_liquid_water_mixing_ratio__dz
         ) / Rq
 
     # pylint: disable=too-many-arguments
     @staticmethod
     def p_of_z_assuming_const_th_and_initial_water_vapour_mixing_ratio(
-        const, g, p0, thstd, water_vapour_mixing_ratio, z
+        const, p0, thstd, water_vapour_mixing_ratio, z
     ):
         z0 = 0
         Rq = const.Rv / (1 / water_vapour_mixing_ratio + 1) + const.Rd / (
             1 + water_vapour_mixing_ratio
         )
         arg = (
             np.power(p0 / const.p1000, const.Rd_over_c_pd)
-            - (z - z0) * const.Rd_over_c_pd * g / thstd / Rq
+            - (z - z0) * const.Rd_over_c_pd * const.g_std / thstd / Rq
         )
         return const.p1000 * np.power(arg, 1 / const.Rd_over_c_pd)

PySDM/physics/hydrostatics/variable_g_isothermal.py

@@ -0,0 +1,25 @@
+"""
+assuming constant temperature and variable
+gravitational acceleration g(z) = g0 * R^2 / (R+z)^2
+as in [Toon et al. 1980](https://doi.org/10.1016/0019-1035(80)90173-6)
+"""
+
+import numpy as np
+from PySDM.physics import constants_defaults
+
+
+class VariableGIsothermal:  # pylint: disable=too-few-public-methods
+    def __init__(self, const):
+        assert np.isfinite(const.celestial_body_radius)
+        assert const.g_std != constants_defaults.g_std
+
+    @staticmethod
+    def pressure(const, z, p0, temperature, molar_mass):
+        return p0 * np.exp(
+            -const.g_std
+            / const.R_str
+            * molar_mass
+            / temperature
+            * z
+            / (1 + z / const.celestial_body_radius)
+        )

PySDM/physics/impl/fake_unit_registry.py

@@ -16,6 +16,7 @@ def __init__(self, si):
             for unit in (
                 "bar",
                 "metre",
+                "meter",
                 "gram",
                 "hertz",
                 "mole",
@@ -25,6 +26,7 @@ def __init__(self, si):
                 "minute",
                 "pascal",
                 "litre",
+                "liter",
                 "hour",
                 "newton",
                 "watt",

docs/bibliography.json

@@ -755,5 +755,13 @@
     ],
     "title": "Isotopic Composition and Temperature of Formation of Antarctic Snows",
     "label": "Picciotto et al. 1960 (Nature 187)"
+  },
+  "https://doi.org/10.1016/0019-1035(80)90173-6": {
+    "usages": [
+      "examples/PySDM_examples/Toon_et_al_1980/fig_1.ipynb",
+      "PySDM/physics/hydrostatics/variable_g_isothermal.py"
+    ],
+    "title": "A physical model of Titan’s clouds",
+    "label": "Toon et al. 1980 (Icarus 43)"
   }
 }

examples/PySDM_examples/Morrison_and_Grabowski_2007/cumulus.py

@@ -50,7 +50,7 @@ def __rhod_of_z(self, T_of_z, p_of_z, q_of_z, z_of_p, z_points):
         def drhod_dz(z, _):
             lv = self.formulae.latent_heat.lv(T_of_z(z))
             return self.formulae.hydrostatics.drho_dz(
-                self.formulae.constants.g_std, p_of_z(z), T_of_z(z), q_of_z(z), lv
+                p_of_z(z), T_of_z(z), q_of_z(z), lv
             )
 
         theta_std0 = self.formulae.trivia.th_std(