shift from tolerance for thd to tolerance for RH in the diffusional growth adaptivity scheme

PySDM/backends/impl_numba/methods/condensation_methods.py

@@ -29,14 +29,15 @@ def condensation(
         rhod,
         thd,
         qv,
+        RH,
         dv,
         prhod,
         pthd,
         pqv,
         kappa,
         f_org,
         rtol_x,
-        rtol_thd,
+        rtol_RH,
         timestep,
         counters,
         cell_order,
@@ -58,14 +59,15 @@ def condensation(
             rhod=rhod.data,
             thd=thd.data,
             qv=qv.data,
+            RH=RH.data,
             dv_mean=dv,
             prhod=prhod.data,
             pthd=pthd.data,
             pqv=pqv.data,
             kappa=kappa.data,
             f_org=f_org.data,
             rtol_x=rtol_x,
-            rtol_thd=rtol_thd,
+            rtol_RH=rtol_RH,
             timestep=timestep,
             counter_n_substeps=counters["n_substeps"].data,
             counter_n_activating=counters["n_activating"].data,
@@ -92,14 +94,15 @@ def _condensation(
         rhod,
         thd,
         qv,
+        RH,
         dv_mean,
         prhod,
         pthd,
         pqv,
         kappa,
         f_org,
         rtol_x,
-        rtol_thd,
+        rtol_RH,
         timestep,
         counter_n_substeps,
         counter_n_activating,
@@ -145,12 +148,13 @@ def _condensation(
                     thd[cell_id],
                     qv[cell_id],
                     rhod[cell_id],
+                    RH[cell_id],
                     dthd_dt,
                     dqv_dt,
                     drhod_dt,
                     md,
                     rtol_x,
-                    rtol_thd,
+                    rtol_RH,
                     timestep,
                     counter_n_substeps[cell_id],
                 )
@@ -177,7 +181,7 @@ def make_adapt_substeps(
         n_substeps_min = math.ceil(timestep / dt_range[1])
 
         @numba.njit(**jit_flags)
-        def adapt_substeps(args, n_substeps, thd, rtol_thd):
+        def adapt_substeps(args, n_substeps, thd, RH, rtol_RH):
             n_substeps = np.maximum(n_substeps_min, n_substeps // multiplier)
             success = False
             for burnout in range(fuse + 1):
@@ -191,24 +195,26 @@ def adapt_substeps(args, n_substeps, thd, rtol_thd):
                         ),
                         return_value=(0, False),
                     )
-                thd_new_long, success = step_fake(args, timestep, n_substeps)
+                _, RH_new_long, success = step_fake(args, timestep, n_substeps)
                 if success:
                     break
                 n_substeps *= multiplier
             for burnout in range(fuse + 1):
                 if burnout == fuse:
                     return warn("burnout (short)", __file__, return_value=(0, False))
-                thd_new_short, success = step_fake(
+                _, RH_new_short, success = step_fake(
                     args, timestep, n_substeps * multiplier
                 )
                 if not success:
                     return warn("short failed", __file__, return_value=(0, False))
-                dthd_long = thd_new_long - thd
-                dthd_short = thd_new_short - thd
-                error_estimate = np.abs(dthd_long - multiplier * dthd_short)
-                thd_new_long = thd_new_short
-                if within_tolerance(error_estimate, thd, rtol_thd):
+
+                dRH_long = RH_new_long - RH
+                dRH_short = RH_new_short - RH
+                error_estimate = np.abs(dRH_long - multiplier * dRH_short)
+                RH_new_long = RH_new_short
+                if within_tolerance(error_estimate, RH, rtol_RH):
                     break
+
                 n_substeps *= multiplier
                 if n_substeps > n_substeps_max:
                     break
@@ -221,8 +227,8 @@ def make_step_fake(jit_flags, step_impl):
         @numba.njit(**jit_flags)
         def step_fake(args, dt, n_substeps):
             dt /= n_substeps
-            _, thd_new, _, _, _, _, success = step_impl(*args, dt, 1, True)
-            return thd_new, success
+            _, thd_new, RH_new, _, _, _, _, success = step_impl(*args, dt, 1, True)
+            return thd_new, RH_new, success
 
         return step_fake
 
@@ -333,6 +339,7 @@ def step_impl(  # pylint: disable=too-many-arguments,too-many-locals
             return (
                 qv,
                 thd,
+                RH,
                 count_activating,
                 count_deactivating,
                 count_ripening,
@@ -417,7 +424,7 @@ def calculate_ml_new(  # pylint: disable=too-many-arguments,too-many-statements,
             lambdaK = phys_lambdaK(T, p)
             lambdaD = phys_lambdaD(DTp, T)
             for drop in cell_idx:
-                if v[drop] < 0:
+                if v[drop] < 0:  # TODO #1086: use _unfrozen from freezing_methods.py
                     continue
                 x_old = x(v[drop])
                 r_old = radius(v[drop])
@@ -666,12 +673,13 @@ def solve(  # pylint: disable=too-many-arguments
             thd,
             qv,
             rhod,
+            RH,
             dthd_dt,
             dqv_dt,
             drhod_dt,
             m_d,
             rtol_x,
-            rtol_thd,
+            rtol_RH,
             timestep,
             n_substeps,
         ):
@@ -694,11 +702,12 @@ def solve(  # pylint: disable=too-many-arguments
             )
             success = True
             if adaptive:
-                n_substeps, success = adapt_substeps(args, n_substeps, thd, rtol_thd)
+                n_substeps, success = adapt_substeps(args, n_substeps, thd, RH, rtol_RH)
             if success:
                 (
                     qv,
                     thd,
+                    _,
                     n_activating,
                     n_deactivating,
                     n_ripening,

PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py

@@ -24,7 +24,7 @@ def patch_particulator(particulator):
 
 
 def _condensation(
-    particulator, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order
+    particulator, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order
 ):
     func = Numba._condensation
     if not numba.config.DISABLE_JIT:  # pylint: disable=no-member
@@ -42,14 +42,15 @@ def _condensation(
         rhod=particulator.environment["rhod"].data,
         thd=particulator.environment["thd"].data,
         qv=particulator.environment["qv"].data,
+        RH=particulator.environment["RH"].data,
         dv_mean=particulator.environment.dv,
         prhod=particulator.environment.get_predicted("rhod").data,
         pthd=particulator.environment.get_predicted("thd").data,
         pqv=particulator.environment.get_predicted("qv").data,
         kappa=particulator.attributes["kappa"].data,
         f_org=particulator.attributes["dry volume organic fraction"].data,
         rtol_x=rtol_x,
-        rtol_thd=rtol_thd,
+        rtol_RH=rtol_RH,
         timestep=particulator.dt,
         counter_n_substeps=counters["n_substeps"],
         counter_n_activating=counters["n_activating"],
@@ -178,14 +179,15 @@ def solve(  # pylint: disable=too-many-arguments,too-many-locals
         thd,
         qv,
         rhod,
+        __,
         dthd_dt,
         dqv_dt,
         drhod_dt,
         m_d_mean,
-        __,
         ___,
-        dt,
         ____,
+        dt,
+        _____,
     ):
         n_sd_in_cell = len(cell_idx)
         y0 = np.empty(n_sd_in_cell + idx_x)

PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py

@@ -332,14 +332,15 @@ def condensation(
         rhod,
         thd,
         qv,
+        RH,
         dv,
         prhod,
         pthd,
         pqv,
         kappa,
         f_org,
         rtol_x,
-        rtol_thd,
+        rtol_RH,
         timestep,
         counters,
         cell_order,

PySDM/dynamics/condensation.py

@@ -8,9 +8,9 @@
 
 from ..physics import si
 
-DEFAULTS = namedtuple("_", ("rtol_x", "rtol_thd", "cond_range", "schedule"))(
+DEFAULTS = namedtuple("_", ("rtol_x", "rtol_RH", "cond_range", "schedule"))(
     rtol_x=1e-6,
-    rtol_thd=1e-6,
+    rtol_RH=1e-6,
     cond_range=(1e-4 * si.second, 1 * si.second),
     schedule="dynamic",
 )
@@ -21,7 +21,7 @@ def __init__(
         self,
         *,
         rtol_x=DEFAULTS.rtol_x,
-        rtol_thd=DEFAULTS.rtol_thd,
+        rtol_RH=DEFAULTS.rtol_RH,
         substeps: int = 1,
         adaptive: bool = True,
         dt_cond_range: tuple = DEFAULTS.cond_range,
@@ -38,7 +38,7 @@ def __init__(
         self.enable = True
 
         self.rtol_x = rtol_x
-        self.rtol_thd = rtol_thd
+        self.rtol_RH = rtol_RH
 
         self.rh_max = None
         self.success = None
@@ -62,7 +62,7 @@ def register(self, builder):
             adaptive=self.adaptive,
             fuse=32,
             multiplier=2,
-            RH_rtol=1e-7,
+            RH_rtol=1e-7,  # TODO #868: use rtol_RH?
             max_iters=self.max_iters,
         )
         builder.request_attribute("critical volume")
@@ -99,7 +99,7 @@ def __call__(self):
 
             self.particulator.condensation(
                 rtol_x=self.rtol_x,
-                rtol_thd=self.rtol_thd,
+                rtol_RH=self.rtol_RH,
                 counters=self.counters,
                 RH_max=self.rh_max,
                 success=self.success,

PySDM/particulator.py

@@ -101,7 +101,7 @@ def update_TpRH(self):
             RH=self.environment.get_predicted("RH"),
         )
 
-    def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order):
+    def condensation(self, *, rtol_x, rtol_RH, counters, RH_max, success, cell_order):
         """Updates droplet volumes by simulating condensation driven by prior changes
           in environment thermodynamic state, updates the environment state.
         In the case of parcel environment, condensation is driven solely by changes in
@@ -121,14 +121,15 @@ def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_orde
             rhod=self.environment["rhod"],
             thd=self.environment["thd"],
             qv=self.environment["qv"],
+            RH=self.environment["RH"],
             dv=self.environment.dv,
             prhod=self.environment.get_predicted("rhod"),
             pthd=self.environment.get_predicted("thd"),
             pqv=self.environment.get_predicted("qv"),
             kappa=self.attributes["kappa"],
             f_org=self.attributes["dry volume organic fraction"],
             rtol_x=rtol_x,
-            rtol_thd=rtol_thd,
+            rtol_RH=rtol_RH,
             v_cr=self.attributes["critical volume"],
             timestep=self.dt,
             counters=counters,

examples/PySDM_examples/Arabas_and_Shima_2017/settings.py

@@ -37,7 +37,7 @@ def __init__(
         self.n_output = 500
 
         self.rtol_x = condensation.DEFAULTS.rtol_x
-        self.rtol_thd = condensation.DEFAULTS.rtol_thd
+        self.rtol_RH = condensation.DEFAULTS.rtol_RH
         self.coord = "volume logarithm"
         self.dt_cond_range = condensation.DEFAULTS.cond_range
 

examples/PySDM_examples/Arabas_and_Shima_2017/simulation.py

@@ -34,7 +34,7 @@ def __init__(self, settings, backend=CPU):
         builder.add_dynamic(
             Condensation(
                 rtol_x=settings.rtol_x,
-                rtol_thd=settings.rtol_thd,
+                rtol_RH=settings.rtol_RH,
                 dt_cond_range=settings.dt_cond_range,
             )
         )

examples/PySDM_examples/Bartman_2020_MasterThesis/fig_5_SCIPY_VS_ADAPTIVE.py

@@ -34,7 +34,7 @@ def data(n_output, rtols, schemes, setups_num):
                 for settings_idx in range(setups_num):
                     settings = setups[settings_idx]
                     settings.rtol_x = rtol
-                    settings.rtol_thd = rtol
+                    settings.rtol_RH = rtol
                     settings.n_output = n_output
                     simulation = Simulation(
                         settings, backend=CPU if scheme == "CPU" else GPU

examples/PySDM_examples/Bartman_et_al_2021/demo.ipynb

@@ -60,7 +60,7 @@
     "settings.simulation_time = .175 * settings.spin_up_time\n",
     "settings.output_interval = 1 * si.minute\n",
     "settings.condensation_rtol_x = 1e-6\n",
-    "settings.condensation_rtol_thd = 5e-7\n",
+    "settings.condensation_rtol_RH = 5e-7\n",
     "\n",
     "settings.condensation_dt_cond_range = (.25*si.s, settings.dt)\n",
     "settings.coalescence_dt_coal_range = settings.condensation_dt_cond_range\n",

examples/PySDM_examples/Bartman_et_al_2021/demo_fig2.ipynb

@@ -45,10 +45,10 @@
     "tol = .5e-6\n",
     "m = 2\n",
     "runs = (\n",
-    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol}},\n",
-    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m}},\n",
-    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m*m}},\n",
-    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_thd': tol*m*m*m}},\n",
+    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol}},\n",
+    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m}},\n",
+    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m*m}},\n",
+    "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': True, 'condensation_rtol_RH': tol*m*m*m}},\n",
     "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps': 128}},\n",
     "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps':  32}},\n",
     "    {'file': TemporaryFile('.nc'), 'settings': {'condensation_adaptive': False, 'condensation_substeps':   8}},\n",
@@ -87,7 +87,7 @@
     "    settings.simulation_time = settings.spin_up_time * (1 if 'CI' not in os.environ else .1)\n",
     "    settings.output_interval = settings.dt\n",
     "    settings.condensation_rtol_x = 1e-6\n",
-    "    settings.condensation_rtol_thd = -1\n",
+    "    settings.condensation_rtol_RH = -1\n",
     "    settings.condensation_schedule = 'dynamic'\n",
     "    settings.kappa = .8\n",
     "    \n",

examples/PySDM_examples/Bartman_et_al_2021/demo_fig3.ipynb

@@ -88,7 +88,7 @@
     "    settings.output_interval = settings.dt\n",
     "    settings.condensation_adaptive = True\n",
     "    settings.condensation_rtol_x = 1e-6\n",
-    "    settings.condensation_rtol_thd = 2e-5/7/7\n",
+    "    settings.condensation_rtol_RH = 2e-5/7/7\n",
     "    settings.condensation_schedule = 'dynamic'\n",
     "    settings.kappa = .8\n",
     "    \n",

examples/PySDM_examples/Morrison_and_Grabowski_2007/common.py

@@ -21,7 +21,7 @@ def __init__(self, formulae: Formulae):
         const = formulae.constants
 
         self.condensation_rtol_x = condensation.DEFAULTS.rtol_x
-        self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd
+        self.condensation_rtol_RH = condensation.DEFAULTS.rtol_RH
         self.condensation_adaptive = True
         self.condensation_substeps = -1
         self.condensation_dt_cond_range = condensation.DEFAULTS.cond_range

examples/PySDM_examples/Pyrcel/simulation.py

@@ -13,7 +13,7 @@
 
 class Simulation(BasicSimulation):
     def __init__(
-        self, settings, products=None, scipy_solver=False, rtol_thd=1e-10, rtol_x=1e-10
+        self, settings, products=None, scipy_solver=False, rtol_RH=1e-6, rtol_x=1e-10
     ):
         env = Parcel(
             dt=settings.timestep,
@@ -27,7 +27,7 @@ def __init__(
         builder = Builder(n_sd=n_sd, backend=CPU(formulae=settings.formulae))
         builder.set_environment(env)
         builder.add_dynamic(AmbientThermodynamics())
-        builder.add_dynamic(Condensation(rtol_thd=rtol_thd, rtol_x=rtol_x))
+        builder.add_dynamic(Condensation(rtol_RH=rtol_RH, rtol_x=rtol_x))
 
         volume = env.mass_of_dry_air / settings.initial_air_density
         attributes = {