removal of time-tracking in parcel env. code (in favour of using Particulator's step counter); removal of unneeded fields: formulae and params (#1509)

slayoo · AgnieszkaZaba · web-flow · commit 2f33cbfcb0ae · 2025-02-08T00:30:26.000+01:00
Co-authored-by: Agnieszka Żaba &lt;azaba@agh.edu.pl&gt;
diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py
@@ -30,52 +30,42 @@ def __init__(
         mixed_phase=False,
         variables: Optional[List[str]] = None,
     ):
-        variables = (variables or []) + ["rhod", "z", "t"]
+        variables = (variables or []) + ["rhod", "z"]
         super().__init__(dt, Mesh.mesh_0d(), variables, mixed_phase=mixed_phase)
 
         self.p0 = p0
         self.initial_water_vapour_mixing_ratio = initial_water_vapour_mixing_ratio
         self.T0 = T0
         self.z0 = z0
         self.mass_of_dry_air = mass_of_dry_air
-
         self.w = w if callable(w) else lambda _: w
-
-        self.formulae = None
         self.delta_liquid_water_mixing_ratio = np.nan
-        self.params = None
 
     @property
     def dv(self):
         rhod_mean = (self.get_predicted("rhod")[0] + self["rhod"][0]) / 2
-        return self.formulae.trivia.volume_of_density_mass(
+        return self.particulator.formulae.trivia.volume_of_density_mass(
             rhod_mean, self.mass_of_dry_air
         )
 
     def register(self, builder):
-        self.formulae = builder.particulator.formulae
-        pd0 = self.formulae.trivia.p_d(self.p0, self.initial_water_vapour_mixing_ratio)
-        rhod0 = self.formulae.state_variable_triplet.rhod_of_pd_T(pd0, self.T0)
-        self.mesh.dv = self.formulae.trivia.volume_of_density_mass(
+        formulae = builder.particulator.formulae
+        pd0 = formulae.trivia.p_d(self.p0, self.initial_water_vapour_mixing_ratio)
+        rhod0 = formulae.state_variable_triplet.rhod_of_pd_T(pd0, self.T0)
+        self.mesh.dv = formulae.trivia.volume_of_density_mass(
             rhod0, self.mass_of_dry_air
         )
 
         Moist.register(self, builder)
 
-        params = (
-            self.initial_water_vapour_mixing_ratio,
-            self.formulae.trivia.th_std(pd0, self.T0),
-            rhod0,
-            self.z0,
-            0,
-        )
-        self["water_vapour_mixing_ratio"][:] = params[0]
-        self["thd"][:] = params[1]
-        self["rhod"][:] = params[2]
-        self["z"][:] = params[3]
-        self["t"][:] = params[4]
+        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"][:] = params[0]
+        self._tmp["water_vapour_mixing_ratio"][
+            :
+        ] = self.initial_water_vapour_mixing_ratio
         self.sync_parcel_vars()
         Moist.sync(self)
         self.notify()
@@ -94,7 +84,7 @@ def init_attributes(
             n_in_dv = np.array([n_in_dv])
 
         attributes = {}
-        dry_volume = self.formulae.trivia.volume(radius=r_dry)
+        dry_volume = self.particulator.formulae.trivia.volume(radius=r_dry)
         attributes["kappa times dry volume"] = dry_volume * kappa
         attributes["multiplicity"] = n_in_dv
         r_wet = equilibrate_wet_radii(
@@ -103,43 +93,44 @@ def init_attributes(
             kappa_times_dry_volume=attributes["kappa times dry volume"],
             rtol=rtol,
         )
-        attributes["volume"] = self.formulae.trivia.volume(radius=r_wet)
+        attributes["volume"] = self.particulator.formulae.trivia.volume(radius=r_wet)
         if include_dry_volume_in_attribute:
             attributes["dry volume"] = dry_volume
         return attributes
 
     def advance_parcel_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
+        formulae = self.particulator.formulae
         T = self["T"][0]
         p = self["p"][0]
-        t = self["t"][0]
 
-        dz_dt = self.w(t + dt / 2)  # "mid-point"
+        dz_dt = self.w((self.particulator.n_steps + 1 / 2) * dt)  # "mid-point"
         water_vapour_mixing_ratio = (
             self["water_vapour_mixing_ratio"][0]
             - self.delta_liquid_water_mixing_ratio / 2
         )
 
-        # derivate evaluated at p_old, T_old, mixrat_mid, w_mid
-        drho_dz = self.formulae.hydrostatics.drho_dz(
-            g=self.formulae.constants.g_std,
+        # 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,
-            lv=self.formulae.latent_heat.lv(T),
+            lv=formulae.latent_heat.lv(T),
             d_liquid_water_mixing_ratio__dz=(
                 self.delta_liquid_water_mixing_ratio / dz_dt / dt
             ),
         )
-        drhod_dz = drho_dz
+        drhod_dz = drho_dz  # TODO #407
 
-        self.particulator.backend.explicit_euler(self._tmp["t"], dt, 1)
         self.particulator.backend.explicit_euler(self._tmp["z"], dt, dz_dt)
         self.particulator.backend.explicit_euler(
             self._tmp["rhod"], dt, dz_dt * drhod_dz
         )
 
-        self.mesh.dv = self.formulae.trivia.volume_of_density_mass(
+        self.mesh.dv = formulae.trivia.volume_of_density_mass(
             (self._tmp["rhod"][0] + self["rhod"][0]) / 2, self.mass_of_dry_air
         )
 
diff --git a/examples/PySDM_examples/Rozanski_and_Sonntag_1982/figs_4_5_6.ipynb b/examples/PySDM_examples/Rozanski_and_Sonntag_1982/figs_4_5_6.ipynb
@@ -74,29 +74,29 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "formulae = Formulae(\n",
+    "FORMULAE = Formulae(\n",
     "    isotope_equilibrium_fractionation_factors='HoritaAndWesolowski1994',\n",
     "    isotope_meteoric_water_line='Dansgaard1964',\n",
     "    isotope_ratio_evolution='MerlivatAndJouzel1979',\n",
     ")\n",
-    "BACKEND = CPU(formulae, override_jit_flags={'parallel': False})\n",
-    "const = formulae.constants\n",
-    "trivia = formulae.trivia\n",
+    "BACKEND = CPU(FORMULAE, override_jit_flags={'parallel': False})\n",
+    "CONST = FORMULAE.constants\n",
+    "TRIVIA = FORMULAE.trivia\n",
     "\n",
     "FIG4_CAPTION_PARAMS = {\n",
-    "    'T_init': const.T0 + 25 * si.K,\n",
+    "    'T_init': CONST.T0 + 25 * si.K,\n",
     "    'P_init': 1000 * si.mbar,\n",
-    "    'RH_init': 80 * const.PER_CENT,\n",
+    "    'RH_init': 80 * CONST.PER_CENT,\n",
     "    # note: R0 in the caption, but negative, hence delta\n",
-    "    'delta_2H_init': -74.7 * const.PER_MILLE,\n",
+    "    'delta_2H_init': -74.7 * CONST.PER_MILLE,\n",
     "    # the \"K\" parameter\n",
     "    'isotope_exchange_factor': 1.,\n",
     "    # the \"N_L\" parameter, note: mixing ratio in the caption, but per-volume units, hence density; assuming 1 kg/m3\n",
     "    'autoconversion_mixrat_threshold': 1 * si.g / si.m**3 / (1 * si.kg / si.m**3)\n",
     "}\n",
     "\n",
     "ARBITRARY_PARAMS = {\n",
-    "    'delta_18O_init': formulae.isotope_meteoric_water_line.d18O_of_d2H(\n",
+    "    'delta_18O_init': FORMULAE.isotope_meteoric_water_line.d18O_of_d2H(\n",
     "        delta_2H=FIG4_CAPTION_PARAMS['delta_2H_init']\n",
     "    ),\n",
     "    'N_SUPER_DROPLETS': 1,\n",
@@ -106,10 +106,10 @@
     "}\n",
     "ARBITRARY_PARAMS['PARCEL_CTOR_ARGS'] = {\n",
     "    'p0': FIG4_CAPTION_PARAMS['P_init'],\n",
-    "    'initial_water_vapour_mixing_ratio': const.eps / (  # TODO #1207: use a physics formula \n",
+    "    'initial_water_vapour_mixing_ratio': CONST.eps / (  # TODO #1207: use a physics formula \n",
     "        FIG4_CAPTION_PARAMS['P_init']\n",
     "        / FIG4_CAPTION_PARAMS['RH_init']\n",
-    "        / formulae.saturation_vapour_pressure.pvs_water(FIG4_CAPTION_PARAMS['T_init'])\n",
+    "        / FORMULAE.saturation_vapour_pressure.pvs_water(FIG4_CAPTION_PARAMS['T_init'])\n",
     "        - 1\n",
     "    ),\n",
     "    'T0': FIG4_CAPTION_PARAMS['T_init'],\n",
@@ -162,11 +162,11 @@
     "        alpha_old = {}\n",
     "        dRv__dt = {}\n",
     "        for isotope in self.isotopes:\n",
-    "            alpha_fun = getattr(formulae.isotope_equilibrium_fractionation_factors, f'alpha_l_{isotope}')\n",
+    "            alpha_fun = getattr(self.particulator.formulae.isotope_equilibrium_fractionation_factors, f'alpha_l_{isotope}')\n",
     "            alpha_old[isotope] = alpha_fun(self[\"T\"][0])\n",
     "            alpha_new = alpha_fun(self._tmp[\"T\"][0])\n",
     "            \n",
-    "            dRv__dt[isotope] = self[f'Rv_{isotope}'][0] * self.formulae.isotope_ratio_evolution.d_Rv_over_Rv(\n",
+    "            dRv__dt[isotope] = self[f'Rv_{isotope}'][0] * self.particulator.formulae.isotope_ratio_evolution.d_Rv_over_Rv(\n",
     "                alpha=alpha_old[isotope],\n",
     "                d_alpha=(alpha_new - alpha_old[isotope]) / self.dt,\n",
     "                n_vapour=self['water_vapour_mixing_ratio'][0],\n",
@@ -231,9 +231,9 @@
     "        builder.add_dynamic(dynamics.Condensation())\n",
     "\n",
     "        for isotope in isotopes:\n",
-    "            builder.particulator.environment[f\"Rv_{isotope}\"][:] = formulae.trivia.isotopic_delta_2_ratio(\n",
+    "            builder.particulator.environment[f\"Rv_{isotope}\"][:] = TRIVIA.isotopic_delta_2_ratio(\n",
     "                delta=(FIG4_CAPTION_PARAMS if isotope == \"2H\" else ARBITRARY_PARAMS)[f'delta_{isotope}_init'],\n",
-    "                reference_ratio=getattr(const, f\"VSMOW_R_{isotope}\")\n",
+    "                reference_ratio=getattr(CONST, f\"VSMOW_R_{isotope}\")\n",
     "            )\n",
     "        \n",
     "        super().__init__(particulator = builder.build(\n",
@@ -12223,7 +12223,7 @@
     "PLOT_KEY = 'd'\n",
     "print(output[PLOT_KEY][0]['Δz'])\n",
     "\n",
-    "temp_C = output[PLOT_KEY][0][\"T\"] - const.T0\n",
+    "temp_C = output[PLOT_KEY][0][\"T\"] - CONST.T0\n",
     "alt_km = in_unit(output[PLOT_KEY][0][\"z\"], si.km) \n",
     "\n",
     "fig, axs = pyplot.subplots(2, 4, squeeze=False, sharey=True, figsize=(11.5, 11.5), tight_layout=True)\n",
@@ -12238,7 +12238,7 @@
     "axs[0,0].yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
     "\n",
     "axs[0,0].plot(\n",
-    "    output[PLOT_KEY][0]['RH'] / const.PER_CENT,\n",
+    "    output[PLOT_KEY][0]['RH'] / CONST.PER_CENT,\n",
     "    temp_C,\n",
     "    color='purple',\n",
     "    label='RH',\n",
@@ -12277,9 +12277,9 @@
     "        (ARBITRARY_PARAMS['N_ITERATIONS']-1) // ARBITRARY_PARAMS['N_PLOT_PROFILES']\n",
     "    ):\n",
     "        deltas = {\n",
-    "            iso: formulae.trivia.isotopic_ratio_2_delta(\n",
+    "            iso: TRIVIA.isotopic_ratio_2_delta(\n",
     "                ratio=output[PLOT_KEY][iteration][f'{ratio}_{iso}'][level_indices['CB']:],\n",
-    "                reference_ratio=getattr(const, f\"VSMOW_R_{iso}\")\n",
+    "                reference_ratio=getattr(CONST, f\"VSMOW_R_{iso}\")\n",
     "            ) for iso in ISOTOPES\n",
     "        }\n",
     "    \n",
@@ -12292,15 +12292,15 @@
     "        for ax, iso, color in ((axs[plot_row, 1], '2H', 'green'), (axs[plot_row, 2], \"18O\", 'brown')):\n",
     "            ax.set_title(title)\n",
     "            ax.plot(\n",
-    "                in_unit(deltas[iso], const.PER_MILLE),\n",
+    "                in_unit(deltas[iso], CONST.PER_MILLE),\n",
     "                temp_C[level_indices['CB']:],\n",
     "                color=color,\n",
     "                **iter_kwargs\n",
     "            )\n",
     "        axs[plot_row, 3].plot(\n",
     "            in_unit(\n",
-    "                formulae.isotope_meteoric_water_line.excess_d(delta_2H=deltas[\"2H\"], delta_18O=deltas[\"18O\"]),\n",
-    "                const.PER_MILLE\n",
+    "                FORMULAE.isotope_meteoric_water_line.excess_d(delta_2H=deltas[\"2H\"], delta_18O=deltas[\"18O\"]),\n",
+    "                CONST.PER_MILLE\n",
     "            ),\n",
     "            temp_C[level_indices['CB']:],\n",
     "            color='black',\n",
@@ -14719,19 +14719,19 @@
     "    colors = {'vapour': 'red', 'rain': 'blue'}\n",
     "    linewidth = 1 + .5 * PARAMS_KEYS.index(key)\n",
     "    line_v = axs.plot(\n",
-    "        in_unit(formulae.trivia.isotopic_ratio_2_delta(\n",
+    "        in_unit(TRIVIA.isotopic_ratio_2_delta(\n",
     "            ratio=data[key]['CT'][\"Rv_2H\"][1:],\n",
-    "            reference_ratio=const.VSMOW_R_2H\n",
-    "        ), const.PER_MILLE),\n",
+    "            reference_ratio=CONST.VSMOW_R_2H\n",
+    "        ), CONST.PER_MILLE),\n",
     "        label=f\"vapour at the cloud top ({in_unit(levels['CT'], si.km)} km)\",\n",
     "        color=colors['vapour'],\n",
     "        linewidth=linewidth\n",
     "    )\n",
     "\n",
-    "    rain_data = in_unit(formulae.trivia.isotopic_ratio_2_delta(\n",
+    "    rain_data = in_unit(TRIVIA.isotopic_ratio_2_delta(\n",
     "        ratio=data[key]['CB'][\"Rr_2H\"][1:],\n",
-    "        reference_ratio=const.VSMOW_R_2H\n",
-    "    ), const.PER_MILLE)\n",
+    "        reference_ratio=CONST.VSMOW_R_2H\n",
+    "    ), CONST.PER_MILLE)\n",
     "    line_r = twinx.plot(\n",
     "        rain_data,\n",
     "        label=f\"rain at the cloud base ({in_unit(levels['CB'], si.km)} km)\",\n",
diff --git a/examples/PySDM_examples/Yang_et_al_2018/simulation.py b/examples/PySDM_examples/Yang_et_al_2018/simulation.py
@@ -69,6 +69,7 @@ def __init__(self, settings, backend=CPU):
             PySDM_products.ActivatedMeanRadius(
                 name="r_act", count_activated=True, count_unactivated=False
             ),
+            PySDM_products.Time(name="t"),
         ]
 
         attributes = environment.init_attributes(
@@ -88,7 +89,8 @@ def save(self, output):
         volume = _sp.attributes["volume"].to_ndarray()
         output["r"].append(self.formulae.trivia.radius(volume=volume))
         output["S"].append(_sp.environment["RH"][cell_id] - 1)
-        for key in ("water_vapour_mixing_ratio", "T", "z", "t"):
+        output["t"].append(_sp.products["t"].get())
+        for key in ("water_vapour_mixing_ratio", "T", "z"):
             output[key].append(_sp.environment[key][cell_id])
         for key in (
             "dt_cond_max",
diff --git a/tests/smoke_tests/parcel_d/rozanski_and_sonntag_1982/test_figs_4_5_6.py b/tests/smoke_tests/parcel_d/rozanski_and_sonntag_1982/test_figs_4_5_6.py
@@ -24,11 +24,11 @@ class TestFigs456:
     @staticmethod
     def test_fig_5_vapour_asymptote(variables):
         delta_vapour_at_the_cloud_top_per_mille = in_unit(
-            variables["formulae"].trivia.isotopic_ratio_2_delta(
+            variables["FORMULAE"].trivia.isotopic_ratio_2_delta(
                 ratio=variables["data"][variables["PLOT_KEY"]]["CT"]["Rv_2H"][1:],
-                reference_ratio=variables["const"].VSMOW_R_2H,
+                reference_ratio=variables["CONST"].VSMOW_R_2H,
             ),
-            variables["const"].PER_MILLE,
+            variables["CONST"].PER_MILLE,
         )
 
         d_delta = np.diff(delta_vapour_at_the_cloud_top_per_mille)
@@ -45,11 +45,11 @@ def test_fig_5_vapour_asymptote(variables):
     @staticmethod
     def test_fig_5_rain_at_the_cloud_base(variables):
         delta_rain_at_the_cloud_base_per_mille = in_unit(
-            variables["formulae"].trivia.isotopic_ratio_2_delta(
+            variables["FORMULAE"].trivia.isotopic_ratio_2_delta(
                 ratio=variables["data"][variables["PLOT_KEY"]]["CB"]["Rr_2H"][1:],
-                reference_ratio=variables["const"].VSMOW_R_2H,
+                reference_ratio=variables["CONST"].VSMOW_R_2H,
             ),
-            variables["const"].PER_MILLE,
+            variables["CONST"].PER_MILLE,
         )
 
         d_delta = np.diff(delta_rain_at_the_cloud_base_per_mille)
