diff --git a/docs/bibliography.json b/docs/bibliography.json index 69ee7458e5..5ea579922c 100644 --- a/docs/bibliography.json +++ b/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" diff --git a/examples/PySDM_examples/Berry_1967/spectrum_plotter.py b/examples/PySDM_examples/Berry_1967/spectrum_plotter.py index 3b92adc5dd..032611f31e 100644 --- a/examples/PySDM_examples/Berry_1967/spectrum_plotter.py +++ b/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) diff --git a/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py b/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py index 20285f0a8a..f989cec287 100644 --- a/examples/PySDM_examples/Shima_et_al_2009/spectrum_plotter.py +++ b/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)), ) diff --git a/tests/smoke_tests/box/shima_et_al_2009/test_convergence.py b/tests/smoke_tests/box/shima_et_al_2009/test_convergence.py new file mode 100644 index 0000000000..f16ec2d552 --- /dev/null +++ b/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")