|
| 1 | +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring |
| 2 | +import matplotlib |
| 3 | +import numpy as np |
| 4 | +from matplotlib import pyplot |
| 5 | +from PySDM_examples.deJong_Mackay_2022 import Settings0D, run_box_breakup |
| 6 | + |
| 7 | +from PySDM.dynamics.collisions.breakup_fragmentations import Straub2010Nf |
| 8 | +from PySDM.dynamics.collisions.coalescence_efficiencies import Straub2010Ec |
| 9 | +from PySDM.physics import si |
| 10 | + |
| 11 | + |
| 12 | +def test_fig_5(plot=False): |
| 13 | + # arrange |
| 14 | + settings = Settings0D(seed=44) |
| 15 | + steps = [0, 30, 60, 180, 540] |
| 16 | + settings._steps = steps # pylint: disable=protected-access |
| 17 | + settings.n_sd = 2**11 |
| 18 | + settings.warn_overflows = False |
| 19 | + settings.radius_bins_edges = np.logspace( |
| 20 | + np.log10(10 * si.um), np.log10(2e3 * si.um), num=32, endpoint=True |
| 21 | + ) |
| 22 | + settings.coal_eff = Straub2010Ec() |
| 23 | + settings.fragmentation = Straub2010Nf(vmin=settings.X0 * 1e-3, nfmax=10) |
| 24 | + |
| 25 | + # act |
| 26 | + (data_x, data_y, _) = run_box_breakup(settings) |
| 27 | + |
| 28 | + # plot |
| 29 | + cmap = matplotlib.cm.get_cmap("viridis") |
| 30 | + for (j, step) in enumerate(steps): |
| 31 | + if j == 0: |
| 32 | + kwargs = {"color": "k", "linestyle": "--", "label": "initial"} |
| 33 | + else: |
| 34 | + kwargs = { |
| 35 | + "color": cmap(j / len(steps)), |
| 36 | + "linestyle": "-", |
| 37 | + "label": f"t = {step}s", |
| 38 | + } |
| 39 | + pyplot.step(data_x, data_y[j] * settings.rho, **kwargs) |
| 40 | + pyplot.xscale("log") |
| 41 | + pyplot.xlabel("particle radius (um)") |
| 42 | + pyplot.ylabel("dm/dlnr (kg/m$^3$ / unit(ln R)") |
| 43 | + pyplot.legend() |
| 44 | + if plot: |
| 45 | + pyplot.show() |
| 46 | + |
| 47 | + # assert |
| 48 | + peaks_expected = { |
| 49 | + 0: (33, 0.018), |
| 50 | + 30: (92, 0.011), |
| 51 | + 60: (305, 0.012), |
| 52 | + 180: (717, 0.015), |
| 53 | + 540: (717, 0.015), |
| 54 | + } |
| 55 | + |
| 56 | + for (j, step) in enumerate(steps): |
| 57 | + print(step) |
| 58 | + peak = np.argmax(data_y[j]) |
| 59 | + np.testing.assert_approx_equal( |
| 60 | + actual=data_x[peak], desired=peaks_expected[step][0], significant=2 |
| 61 | + ) |
| 62 | + np.testing.assert_approx_equal( |
| 63 | + actual=data_y[j][peak] * settings.rho, |
| 64 | + desired=peaks_expected[step][1], |
| 65 | + significant=2, |
| 66 | + ) |
0 commit comments