refactor saturation vapour pressure formulae to avoid using temperatu…

…res in Celsius (#1388) Co-authored-by: Sylwester Arabas <[email protected]>
open-atmos · Sep 17, 2024 · 80043a8 · 80043a8
1 parent 20bed75
commit 80043a8
Show file tree

Hide file tree

Showing 34 changed files with 3,477 additions and 289 deletions.
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
@@ -18,4 +18,4 @@ repos:
     hooks:
       - id: trailing-whitespace
       - id: end-of-file-fixer
-      - id: debug-statements
+      - id: debug-statements
diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py
@@ -291,9 +291,7 @@ def step_impl(  # pylint: disable=too-many-arguments,too-many-locals
                 )
                 pv = formulae.state_variable_triplet__pv(p, water_vapour_mixing_ratio)
                 lv = formulae.latent_heat__lv(T)
-                pvs = formulae.saturation_vapour_pressure__pvs_Celsius(
-                    T - formulae.constants.T0
-                )
+                pvs = formulae.saturation_vapour_pressure__pvs_water(T)
                 DTp = formulae.diffusion_thermics__D(T, p)
                 RH = pv / pvs
                 Sc = formulae.trivia__air_schmidt_number(

diff --git a/PySDM/backends/impl_numba/methods/physics_methods.py b/PySDM/backends/impl_numba/methods/physics_methods.py
@@ -59,7 +59,7 @@ def body(*, rhod, thd, water_vapour_mixing_ratio, T, p, RH):
                 )
                 RH[i] = ff.state_variable_triplet__pv(
                     p[i], water_vapour_mixing_ratio[i]
-                ) / ff.saturation_vapour_pressure__pvs_Celsius(T[i] - ff.constants.T0)
+                ) / ff.saturation_vapour_pressure__pvs_water(T[i])
 
         return body
 
@@ -82,9 +82,7 @@ def _a_w_ice_body(self):
         @numba.njit(**self.default_jit_flags)
         def body(*, T_in, p_in, RH_in, water_vapour_mixing_ratio_in, a_w_ice_out):
             for i in prange(T_in.shape[0]):  # pylint: disable=not-an-iterable
-                pvi = ff.saturation_vapour_pressure__ice_Celsius(
-                    T_in[i] - ff.constants.T0
-                )
+                pvi = ff.saturation_vapour_pressure__pvs_ice(T_in[i])
                 pv = ff.state_variable_triplet__pv(
                     p_in[i], water_vapour_mixing_ratio_in[i]
                 )

diff --git a/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py b/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py
@@ -18,7 +18,7 @@
     _CellData,
     _Counters,
 )
-from PySDM.physics.constants_defaults import PI_4_3, T0
+from PySDM.physics.constants_defaults import PI_4_3
 
 idx_thd = 0
 idx_x = 1
@@ -205,7 +205,7 @@ def _odesys(  # pylint: disable=too-many-arguments,too-many-locals
         T = jit_formulae.state_variable_triplet__T(rhod, thd)
         p = jit_formulae.state_variable_triplet__p(rhod, T, water_vapour_mixing_ratio)
         pv = jit_formulae.state_variable_triplet__pv(p, water_vapour_mixing_ratio)
-        pvs = jit_formulae.saturation_vapour_pressure__pvs_Celsius(T - T0)
+        pvs = jit_formulae.saturation_vapour_pressure__pvs_water(T)
         RH = pv / pvs
 
         dy_dt = np.empty_like(y)

diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py
@@ -318,8 +318,8 @@ def __pre(self):
                 p='p[i]', water_vapour_mixing_ratio='predicted_water_vapour_mixing_ratio[i]')};
             lv[i] = {phys.latent_heat.lv.c_inline(
                 T='T[i]')};
-            pvs[i] = {phys.saturation_vapour_pressure.pvs_Celsius.c_inline(
-                T='T[i] - const.T0')};
+            pvs[i] = {phys.saturation_vapour_pressure.pvs_water.c_inline(
+                T='T[i]')};
             RH[i] = pv[i] / pvs[i];
             RH_max[i] = max(RH_max[i], RH[i]);
             DTp[i] = {phys.diffusion_thermics.D.c_inline(

diff --git a/PySDM/backends/impl_thrust_rtc/methods/physics_methods.py b/PySDM/backends/impl_thrust_rtc/methods/physics_methods.py
@@ -27,9 +27,7 @@ def _temperature_pressure_rh_body(self):
             )};
             RH[i] = {self.formulae.state_variable_triplet.pv.c_inline(
                 p="p[i]", water_vapour_mixing_ratio="water_vapour_mixing_ratio[i]"
-            )} / {self.formulae.saturation_vapour_pressure.pvs_Celsius.c_inline(
-                T="T[i] - const.T0"
-            )};
+            )} / {self.formulae.saturation_vapour_pressure.pvs_water.c_inline(T="T[i]")};
             """.replace(
                 "real_type", self._get_c_type()
             ),

diff --git a/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py b/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py
@@ -12,10 +12,12 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def pvs_Celsius(const, T):
-        return const.ARM_C1 * np.exp((const.ARM_C2 * T) / (T + const.ARM_C3))
+    def pvs_water(const, T):
+        return const.ARM_C1 * np.exp(
+            (const.ARM_C2 * (T - const.T0)) / ((T - const.T0) + const.ARM_C3)
+        )
 
     @staticmethod
-    def ice_Celsius(const, T):
+    def pvs_ice(const, T):
         """NaN with unit of pressure and correct dimension"""
-        return np.nan * T / const.ARM_C3 * const.ARM_C1
+        return np.nan * (T - const.T0) / const.ARM_C3 * const.ARM_C1
diff --git a/PySDM/physics/saturation_vapour_pressure/bolton_1980.py b/PySDM/physics/saturation_vapour_pressure/bolton_1980.py
@@ -10,11 +10,13 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def pvs_Celsius(const, T):
-        """valid for -30 <= T <= 35 C, eq (10)"""
-        return const.B80W_G0 * np.exp((const.B80W_G1 * T) / (T + const.B80W_G2))
+    def pvs_water(const, T):
+        """valid for 243.15(-30) <= T <= 308.15(35) K(C), eq. (10)"""
+        return const.B80W_G0 * np.exp(
+            (const.B80W_G1 * (T - const.T0)) / ((T - const.T0) + const.B80W_G2)
+        )
 
     @staticmethod
-    def ice_Celsius(const, T):
+    def pvs_ice(const, T):
         """NaN with unit of pressure and correct dimension"""
-        return np.nan * T / const.B80W_G2 * const.B80W_G0
+        return np.nan * (T - const.T0) / const.B80W_G2 * const.B80W_G0
diff --git a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py
@@ -9,45 +9,55 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def pvs_Celsius(const, T):
-        return const.FWC_C0 + T * (
+    def pvs_water(const, T):
+        return const.FWC_C0 + (T - const.T0) * (
             const.FWC_C1
-            + T
+            + (T - const.T0)
             * (
                 const.FWC_C2
-                + T
+                + (T - const.T0)
                 * (
                     const.FWC_C3
-                    + T
+                    + (T - const.T0)
                     * (
                         const.FWC_C4
-                        + T
+                        + (T - const.T0)
                         * (
                             const.FWC_C5
-                            + T * (const.FWC_C6 + T * (const.FWC_C7 + T * const.FWC_C8))
+                            + (T - const.T0)
+                            * (
+                                const.FWC_C6
+                                + (T - const.T0)
+                                * (const.FWC_C7 + (T - const.T0) * const.FWC_C8)
+                            )
                         )
                     )
                 )
             )
         )
 
     @staticmethod
-    def ice_Celsius(const, T):
-        return const.FWC_I0 + T * (
+    def pvs_ice(const, T):
+        return const.FWC_I0 + (T - const.T0) * (
             const.FWC_I1
-            + T
+            + (T - const.T0)
             * (
                 const.FWC_I2
-                + T
+                + (T - const.T0)
                 * (
                     const.FWC_I3
-                    + T
+                    + (T - const.T0)
                     * (
                         const.FWC_I4
-                        + T
+                        + (T - const.T0)
                         * (
                             const.FWC_I5
-                            + T * (const.FWC_I6 + T * (const.FWC_I7 + T * const.FWC_I8))
+                            + (T - const.T0)
+                            * (
+                                const.FWC_I6
+                                + (T - const.T0)
+                                * (const.FWC_I7 + (T - const.T0) * const.FWC_I8)
+                            )
                         )
                     )
                 )

diff --git a/PySDM/physics/saturation_vapour_pressure/lowe1977.py b/PySDM/physics/saturation_vapour_pressure/lowe1977.py
@@ -9,31 +9,41 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def pvs_Celsius(const, T):
-        return const.L77W_A0 + T * (
+    def pvs_water(const, T):
+        return const.L77W_A0 + (T - const.T0) * (
             const.L77W_A1
-            + T
+            + (T - const.T0)
             * (
                 const.L77W_A2
-                + T
+                + (T - const.T0)
                 * (
                     const.L77W_A3
-                    + T * (const.L77W_A4 + T * (const.L77W_A5 + T * (const.L77W_A6)))
+                    + (T - const.T0)
+                    * (
+                        const.L77W_A4
+                        + (T - const.T0)
+                        * (const.L77W_A5 + (T - const.T0) * (const.L77W_A6))
+                    )
                 )
             )
         )
 
     @staticmethod
-    def ice_Celsius(const, T):
-        return const.L77I_A0 + T * (
+    def pvs_ice(const, T):
+        return const.L77I_A0 + (T - const.T0) * (
             const.L77I_A1
-            + T
+            + (T - const.T0)
             * (
                 const.L77I_A2
-                + T
+                + (T - const.T0)
                 * (
                     const.L77I_A3
-                    + T * (const.L77I_A4 + T * (const.L77I_A5 + T * (const.L77I_A6)))
+                    + (T - const.T0)
+                    * (
+                        const.L77I_A4
+                        + (T - const.T0)
+                        * (const.L77I_A5 + (T - const.T0) * (const.L77I_A6))
+                    )
                 )
             )
         )
diff --git a/PySDM/physics/saturation_vapour_pressure/murphy_koop_2005.py b/PySDM/physics/saturation_vapour_pressure/murphy_koop_2005.py
@@ -10,28 +10,28 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def pvs_Celsius(const, T):
+    def pvs_water(const, T):
         """valid for 123 < T < 332 K, eq (10)"""
         return const.MK05_LIQ_C1 * np.exp(
             const.MK05_LIQ_C2
-            - const.MK05_LIQ_C3 / (T + const.T0)
-            - const.MK05_LIQ_C4 * np.log((T + const.T0) / const.MK05_LIQ_C5)
-            + const.MK05_LIQ_C6 * (T + const.T0)
-            + np.tanh(const.MK05_LIQ_C7 * (T + const.T0 - const.MK05_LIQ_C8))
+            - const.MK05_LIQ_C3 / (T)
+            - const.MK05_LIQ_C4 * np.log(T / const.MK05_LIQ_C5)
+            + const.MK05_LIQ_C6 * (T)
+            + np.tanh(const.MK05_LIQ_C7 * (T - const.MK05_LIQ_C8))
             * (
                 const.MK05_LIQ_C9
-                - const.MK05_LIQ_C10 / (T + const.T0)
-                - const.MK05_LIQ_C11 * np.log((T + const.T0) / const.MK05_LIQ_C12)
-                + const.MK05_LIQ_C13 * (T + const.T0)
+                - const.MK05_LIQ_C10 / T
+                - const.MK05_LIQ_C11 * np.log(T / const.MK05_LIQ_C12)
+                + const.MK05_LIQ_C13 * T
             )
         )
 
     @staticmethod
-    def ice_Celsius(const, T):
+    def pvs_ice(const, T):
         """valid for T > 110 K, eq (7)"""
         return const.MK05_ICE_C1 * np.exp(
             const.MK05_ICE_C2
-            - const.MK05_ICE_C3 / (T + const.T0)
-            + const.MK05_ICE_C4 * np.log((T + const.T0) / const.MK05_ICE_C5)
-            - const.MK05_ICE_C6 * (T + const.T0)
+            - const.MK05_ICE_C3 / T
+            + const.MK05_ICE_C4 * np.log(T / const.MK05_ICE_C5)
+            - const.MK05_ICE_C6 * T
         )
diff --git a/PySDM/physics/saturation_vapour_pressure/wexler_1976.py b/PySDM/physics/saturation_vapour_pressure/wexler_1976.py
@@ -10,22 +10,22 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def pvs_Celsius(const, T):
+    def pvs_water(const, T):
         return (
             np.exp(
-                const.W76W_G0 / (T + const.T0) ** 2
-                + const.W76W_G1 / (T + const.T0)
+                const.W76W_G0 / T**2
+                + const.W76W_G1 / T
                 + const.W76W_G2
-                + const.W76W_G3 * (T + const.T0)
-                + const.W76W_G4 * (T + const.T0) ** 2
-                + const.W76W_G5 * (T + const.T0) ** 3
-                + const.W76W_G6 * (T + const.T0) ** 4
-                + const.W76W_G7 * np.log((T + const.T0) / const.one_kelvin)
+                + const.W76W_G3 * T
+                + const.W76W_G4 * T**2
+                + const.W76W_G5 * T**3
+                + const.W76W_G6 * T**4
+                + const.W76W_G7 * np.log(T / const.one_kelvin)
             )
             * const.W76W_G8
         )
 
     @staticmethod
-    def ice_Celsius(const, T):
+    def pvs_ice(const, T):
         """NaN with unit of pressure and correct dimension"""
-        return np.nan * T / const.B80W_G2 * const.B80W_G0
+        return np.nan * (T - const.T0) / const.B80W_G2 * const.B80W_G0
diff --git a/examples/PySDM_examples/Abdul_Razzak_Ghan_2000/run_ARG_parcel.py b/examples/PySDM_examples/Abdul_Razzak_Ghan_2000/run_ARG_parcel.py
@@ -35,7 +35,7 @@ def run_parcel(
 
     formulae = Formulae(constants=CONSTANTS_ARG)
     const = formulae.constants
-    pv0 = RH0 * formulae.saturation_vapour_pressure.pvs_Celsius(T0 - const.T0)
+    pv0 = RH0 * formulae.saturation_vapour_pressure.pvs_water(T0)
 
     env = Parcel(
         dt=dt,

diff --git a/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py b/examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py
@@ -11,7 +11,6 @@
 from PySDM.environments import Box
 from PySDM.initialisation import discretise_multiplicities
 from PySDM.initialisation.sampling import spectral_sampling
-from PySDM.physics import constants as const
 from PySDM.physics import si
 from PySDM.products import IceWaterContent, TotalUnfrozenImmersedSurfaceArea
 
@@ -143,7 +142,7 @@ def plot_j_het(self, variant: str, abifm_params_case: str, ylim=None):
         svp = formulae.saturation_vapour_pressure
         plot_x = np.linspace(*self.temperature_range) * si.K
         plot_y = formulae.heterogeneous_ice_nucleation_rate.j_het(
-            svp.ice_Celsius(plot_x - const.T0) / svp.pvs_Celsius(plot_x - const.T0)
+            svp.pvs_ice(plot_x) / svp.pvs_water(plot_x)
         )
         pyplot.grid()
         pyplot.plot(plot_x, plot_y / yunit, color="red", label="ABIFM $J_{het}$")
@@ -252,9 +251,7 @@ def simulation(
     for i in range(int(total_time / time_step) + 1):
         if cooling_rate != 0:
             env["T"] -= np.full((1,), cooling_rate * time_step / 2)
-            env["a_w_ice"] = svp.ice_Celsius(env["T"][0] - const.T0) / svp.pvs_Celsius(
-                env["T"][0] - const.T0
-            )
+            env["a_w_ice"] = svp.pvs_ice(env["T"][0]) / svp.pvs_water(env["T"][0])
         particulator.run(0 if i == 0 else 1)
         if cooling_rate != 0:
             env["T"] -= np.full((1,), cooling_rate * time_step / 2)

diff --git a/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py b/examples/PySDM_examples/Arabas_and_Shima_2017/settings.py
@@ -27,7 +27,7 @@ def __init__(
         self.T0 = 300 * si.kelvin
         self.z_half = 150 * si.metres
 
-        pvs = self.formulae.saturation_vapour_pressure.pvs_Celsius(self.T0 - const.T0)
+        pvs = self.formulae.saturation_vapour_pressure.pvs_water(self.T0)
         self.initial_water_vapour_mixing_ratio = const.eps / (
             self.p0 / self.RH0 / pvs - 1
         )

diff --git a/examples/PySDM_examples/Arabas_et_al_2023/fig_2.ipynb b/examples/PySDM_examples/Arabas_et_al_2023/fig_2.ipynb
diff --git a/examples/PySDM_examples/Arabas_et_al_2023/run_simulation.py b/examples/PySDM_examples/Arabas_et_al_2023/run_simulation.py
@@ -4,10 +4,9 @@
 def update_thermo(particulator, T):
     env = particulator.environment
     svp = particulator.formulae.saturation_vapour_pressure
-    const = particulator.formulae.constants
 
     env["T"] = T
-    env["a_w_ice"] = svp.ice_Celsius(T - const.T0) / svp.pvs_Celsius(T - const.T0)
+    env["a_w_ice"] = svp.pvs_ice(T) / svp.pvs_water(T)
 
 
 def run_simulation(particulator, temperature_profile, n_steps):

diff --git a/examples/PySDM_examples/Gedzelman_and_Arnold_1994/fig_2.ipynb b/examples/PySDM_examples/Gedzelman_and_Arnold_1994/fig_2.ipynb
@@ -174,7 +174,7 @@
     "        )(kwargs['temperature'])\n",
     "\n",
     "        missing_b_multiplier = (\n",
-    "            kwargs['formulae'].saturation_vapour_pressure.pvs_Celsius(kwargs['temperature'] - T0)\n",
+    "            kwargs['formulae'].saturation_vapour_pressure.pvs_water(kwargs['temperature'])\n",
     "            / kwargs['temperature']\n",
     "            / const.Rv\n",
     "        )\n",

diff --git a/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/settings.py b/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/settings.py
@@ -61,9 +61,7 @@ def __init__(
         self.initial_temperature = initial_temperature
         pv0 = (
             initial_relative_humidity
-            * self.formulae.saturation_vapour_pressure.pvs_Celsius(
-                initial_temperature - const.T0
-            )
+            * self.formulae.saturation_vapour_pressure.pvs_water(initial_temperature)
         )
         self.initial_vapour_mixing_ratio = const.eps * pv0 / (initial_pressure - pv0)
         self.t_max = displacement / vertical_velocity