add smoke test checking error measures for the Shima_et_al_2009 box example. Closes #327 (#1191)

Co-authored-by: AgnieszkaZaba <56157996+AgnieszkaZaba@users.noreply.github.com>
Co-authored-by: Agnieszka Żaba <azaba@agh.edu.pl>

Author: slayoo

docs/bibliography.json

@@ -121,7 +121,8 @@
       "examples/PySDM_examples/Shima_et_al_2009/__init__.py",
       "PySDM/dynamics/collisions/collision.py",
       "tutorials/collisions/collisions_playground.ipynb",
-      "tutorials/wikipedia/sdm.ipynb"
+      "tutorials/wikipedia/sdm.ipynb",
+      "tests/smoke_tests/box/shima_et_al_2009/test_convergence.py"
     ],
     "label": "Shima et al. 2009 (Q. J. R. Meteorol. Soc. 135)",
     "title": "The super-droplet method for the numerical simulation of clouds and precipitation: a particle-based and probabilistic microphysics model coupled with a non-hydrostatic model"

examples/PySDM_examples/Berry_1967/spectrum_plotter.py

@@ -26,6 +26,8 @@ def show(self):
         self.ticks()
         show_plot()
 
-    def plot(self, spectrum, t):
+    def plot(
+        self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False
+    ):
         settings = self.settings
-        self.plot_data(settings, t, spectrum)
+        self.plot_data(settings, t, spectrum, label, color)

examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py

@@ -81,12 +81,16 @@ def save(self, file):
         self.finish()
         pyplot.savefig(file, format=self.format)
 
-    def plot(self, spectrum, t):
-        error = self.plot_analytic_solution(self.settings, t, spectrum)
-        self.plot_data(self.settings, t, spectrum)
+    def plot(
+        self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False
+    ):
+        error = self.plot_analytic_solution(self.settings, t, spectrum, title)
+        if label is not None and add_error_to_label:
+            label += f" error={error:.4g}"
+        self.plot_data(self.settings, t, spectrum, label, color)
         return error
 
-    def plot_analytic_solution(self, settings, t, spectrum=None):
+    def plot_analytic_solution(self, settings, t, spectrum, title):
         if t == 0:
             analytic_solution = settings.spectrum.size_distribution
         else:
@@ -123,11 +127,13 @@ 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 = (
+                title or f"error measure: {error:.2f}"
+            )  # TODO #327 relative error
             return error
         return None
 
-    def plot_data(self, settings, t, spectrum):
+    def plot_data(self, settings, t, spectrum, label, color):
         if self.smooth:
             scope = self.smooth_scope
             if t != 0:
@@ -144,18 +150,16 @@ def plot_data(self, settings, t, spectrum):
             self.ax.plot(
                 (x[:-1] + dx / 2) * si.metres / si.micrometres,
                 spectrum[:-scope] * si.kilograms / si.grams,
-                label=f"t = {t}s",
-                color=self.colors(
-                    t / (self.settings.output_steps[-1] * self.settings.dt)
-                ),
+                label=label or f"t = {t}s",
+                color=color
+                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
             )
         else:
             self.ax.step(
                 settings.radius_bins_edges[:-1] * si.metres / si.micrometres,
                 spectrum * si.kilograms / si.grams,
                 where="post",
-                label=f"t = {t}s",
-                color=self.colors(
-                    t / (self.settings.output_steps[-1] * self.settings.dt)
-                ),
+                label=label or f"t = {t}s",
+                color=color
+                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
             )

tests/smoke_tests/box/shima_et_al_2009/test_convergence.py

@@ -0,0 +1,85 @@
+"""
+Tests checking if, for the example case from Fig 4 in
+[Shima et al. 2009](https://doi.org/10.1002/qj.441),
+the simulations converge towards analytic solution.
+"""
+
+from inspect import signature
+from itertools import islice
+
+import pytest
+from matplotlib import pyplot
+from PySDM_examples.Shima_et_al_2009.example import run
+from PySDM_examples.Shima_et_al_2009.settings import Settings
+from PySDM_examples.Shima_et_al_2009.spectrum_plotter import SpectrumPlotter
+
+from PySDM.physics import si
+
+COLORS = ("red", "green", "blue")
+
+
+class TestConvergence:  # pylint: disable=missing-class-docstring
+    @staticmethod
+    @pytest.mark.parametrize(
+        "adaptive, dt",
+        (
+            pytest.param(False, 100 * si.s, marks=pytest.mark.xfail(strict=True)),
+            (True, 100 * si.s),
+            pytest.param(False, 50 * si.s, marks=pytest.mark.xfail(strict=True)),
+            (True, 50 * si.s),
+        ),
+    )
+    def test_convergence_with_sd_count(dt, adaptive, plot=False):
+        """check if increasing the number of super particles indeed
+        reduces the error of the simulation (vs. analytic solution)"""
+        # arrange
+        settings = Settings(steps=[3600])
+        settings.adaptive = adaptive
+        plotter = SpectrumPlotter(settings)
+        errors = {}
+
+        # act
+        for i, ln2_nsd in enumerate((11, 15, 19)):
+            settings.dt = dt
+            settings.n_sd = 2**ln2_nsd
+            values, _ = run(settings)
+
+            title = (
+                ""
+                if i != 0
+                else (
+                    f"{settings.dt=}  settings.times={settings.steps}  {settings.adaptive=}"
+                )
+            )
+            errors[ln2_nsd] = plotter.plot(
+                **dict(
+                    islice(
+                        {  # supporting older versions of PySDM-examples
+                            "t": settings.steps[-1],
+                            "spectrum": values[tuple(values.keys())[-1]],
+                            "label": f"{ln2_nsd=}",
+                            "color": COLORS[i],
+                            "title": title,
+                            "add_error_to_label": True,
+                        }.items(),
+                        len(signature(plotter.plot).parameters),
+                    )
+                )
+            )
+
+        # plot
+        if plot:
+            plotter.show()
+        else:
+            # https://github.com/matplotlib/matplotlib/issues/9970
+            pyplot.xscale("linear")
+            pyplot.clf()
+
+        # assert monotonicity (i.e., the larger the sd count, the smaller the error)
+        assert tuple(errors.keys()) == tuple(sorted(errors.keys()))
+        assert tuple(errors.values()) == tuple(reversed(sorted(errors.values())))
+
+    @staticmethod
+    def test_convergence_with_timestep():
+        """ditto for timestep"""
+        pytest.skip("# TODO #1189")