adaptivity examples beef up

PySDM/backends/impl_common/backend_methods.py

@@ -3,8 +3,11 @@
  class for all backend methods classes
 """
 
+from pystrict import strict
+
 
 # pylint: disable=too-few-public-methods
+@strict
 class BackendMethods:
     def __init__(self):
         if not hasattr(self, "formulae"):

PySDM/backends/impl_numba/methods/physics_methods.py

@@ -125,6 +125,17 @@
 
         return body
 
+    @cached_property
+    def mass_of_volume(self):
+        phys_volume_to_mass = self.formulae.particle_shape_and_density.volume_to_mass
+
+        @numba.njit(**self.default_jit_flags)
+        def mass_of_volume(mass, volume):
+            for i in prange(volume.shape[0]):  # pylint: disable=not-an-iterable
+                mass[i] = phys_volume_to_mass(volume[i])
+
+        return mass_of_volume
+
     def mass_of_water_volume(self, mass, volume):
         self._mass_of_volume_body(mass.data, volume.data)
 

examples/PySDM_examples/Bartman_2020_MasterThesis/fig_4_adaptive_sdm.py

@@ -1,4 +1,5 @@
 import os
+import numpy as np
 
 from matplotlib import pyplot as plt
 from PySDM_examples.Shima_et_al_2009.example import run
@@ -7,10 +8,11 @@
 
 
 def main(plot: bool = True, save: str = None):
-    n_sds = [13, 15, 17]
-    dts = [10, 5, 1, "adaptive"]
-    iters = 10
+    n_sds = [13, 15, 17, 19]
+    dts = [10, 1, "adaptive"]
+    iters_without_warmup = 5
     base_time = None
+    base_error = None
 
     plt.ioff()
     fig, axs = plt.subplots(
@@ -20,17 +22,20 @@ def main(plot: bool = True, save: str = None):
     for i, dt in enumerate(dts):
         for j, n_sd in enumerate(n_sds):
             outputs = []
-            exec_time = 0
-            for _ in range(iters):
-                settings = Settings()
+            exec_times = []
+            one_for_warmup = 1
+            for it in range(iters_without_warmup + one_for_warmup):
+                settings = Settings(seed=it)
 
                 settings.n_sd = 2**n_sd
-                settings.dt = dt if dt != "adaptive" else 10
+                settings.dt = dt if dt != "adaptive" else max(dts[:-1])
                 settings.adaptive = dt == "adaptive"
 
                 states, exec_time = run(settings)
+                print(f"{dt=}, {n_sd=}, {exec_time=}, {it=}")
+                exec_times.append(exec_time)
                 outputs.append(states)
-            mean_time = exec_time / iters
+            mean_time = np.mean(exec_times[one_for_warmup:])
             if base_time is None:
                 base_time = mean_time
             norm_time = mean_time / base_time
@@ -45,10 +50,17 @@ def main(plot: bool = True, save: str = None):
             plotter.ax = axs[i, j]
             plotter.smooth = True
             for step, vals in mean_output.items():
-                plotter.plot(vals, step * settings.dt)
+                error = plotter.plot(vals, step * settings.dt)
+            last_step_error = error
+            if base_error is None:
+                base_error = last_step_error
+            norm_error = last_step_error / base_error
 
             plotter.ylabel = (
-                r"$\bf{dt: " + str(dt) + "}$\ndm/dlnr [g/m^3/(unit dr/r)]"
+                r"$\bf{dt: "
+                + str(settings.dt)
+                + ("+ adaptivity" if settings.adaptive else "")
+                + "}$\ndm/dlnr [g/m^3/(unit dr/r)]"
                 if j == 0
                 else None
             )
@@ -57,7 +69,9 @@ def main(plot: bool = True, save: str = None):
                 if i == len(dts) - 1
                 else None
             )
-            plotter.title = f"norm. time: {norm_time:.2f}; " + plotter.title
+            plotter.title = (
+                f"time: {norm_time:.2f} error: {norm_error:.2f} (normalised)"
+            )
             plotter.finished = False
             plotter.finish()
     if save is not None:

examples/PySDM_examples/Shima_et_al_2009/error_measure.py

@@ -1,11 +1,15 @@
 import numpy as np
 
 
-def error_measure(y, y_true, x):
+def error_measure_old(y, y_true, x):
     errors = y_true - y
     errors = errors[0:-1] + errors[1:]
     dx = np.diff(x)
     errors *= dx
     errors /= 2
     error = np.sum(np.abs(errors))
     return error
+
+
+def error_measure(y, y_true, _):
+    return np.sqrt(np.mean(np.square(y - y_true)))

examples/PySDM_examples/Shima_et_al_2009/example.py

@@ -10,7 +10,7 @@
 from PySDM.dynamics import Coalescence
 from PySDM.environments import Box
 from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
-from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum, WallTime
+from PySDM.products import ParticleVolumeVersusRadiusLogarithmSpectrum, CPUTime
 
 
 def run(settings, backend=CPU, observers=()):
@@ -29,21 +29,21 @@ def run(settings, backend=CPU, observers=()):
         ParticleVolumeVersusRadiusLogarithmSpectrum(
             settings.radius_bins_edges, name="dv/dlnr"
         ),
-        WallTime(),
+        CPUTime(),
     )
     particulator = builder.build(attributes, products)
 
     for observer in observers:
         particulator.observers.append(observer)
 
     vals = {}
-    particulator.products["wall time"].reset()
+    particulator.products["CPU Time"].reset()
     for step in settings.output_steps:
         particulator.run(step - particulator.n_steps)
         vals[step] = particulator.products["dv/dlnr"].get()[0]
         vals[step][:] *= settings.rho
 
-    exec_time = particulator.products["wall time"].get()
+    exec_time = particulator.products["CPU Time"].get()
     return vals, exec_time
 
 

examples/PySDM_examples/Shima_et_al_2009/settings.py

@@ -11,9 +11,9 @@
 
 @strict
 class Settings:
-    def __init__(self, steps: Optional[list] = None):
+    def __init__(self, seed: Optional[int] = None, steps: Optional[list] = None):
         steps = steps or [0, 1200, 2400, 3600]
-        self.formulae = Formulae()
+        self.formulae = Formulae(seed=seed)
         self.n_sd = 2**13
         self.n_part = 2**23 / si.metre**3
         self.X0 = self.formulae.trivia.volume(radius=30.531 * si.micrometres)
@@ -22,7 +22,6 @@ def __init__(self, steps: Optional[list] = None):
         self.rho = 1000 * si.kilogram / si.metre**3
         self.dt = 1 * si.seconds
         self.adaptive = False
-        self.seed = 44
         self.steps = steps
         self.kernel = Golovin(b=1.5e3 / si.second)
         self.spectrum = spectra.Exponential(norm_factor=self.norm_factor, scale=self.X0)

examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py

@@ -123,7 +123,8 @@ def analytic_solution(x):
         if spectrum is not None:
             y = spectrum * si.kilograms / si.grams
             error = error_measure(y, y_true, x)
-            self.title = f"error measure: {error:.2f}"  # TODO #327 relative error
+            self.title = f"error: {error:.2f}"  # TODO #327 relative error
+            self.title = f"error: {error:.2f}"  # TODO #327 relative error
             return error
         return None
 

setup.py

@@ -69,6 +69,7 @@ def get_long_description():
             else ""
         ),
         "pyevtk" + ("==1.2.0" if CI else ""),
+        "pystrict",
     ],
     extras_require={
         "tests": [