Merge pull request #616 from slayoo/freezing

adding smoke test for Alpert & Knopf 2016 Fig 1

Author: slayoo

PySDM/physics/constants.py

@@ -66,6 +66,16 @@ def convert_to(value, unit):
 FWC_C7 = 0.379534310e-13 * si.hPa / si.K**7
 FWC_C8 = -.321582393e-15 * si.hPa / si.K**8
 
+FWC_I0 = 6.098689930e000 * si.hPa
+FWC_I1 = 0.499320233e000 * si.hPa / si.K
+FWC_I2 = 0.184672631e-01 * si.hPa / si.K**2
+FWC_I3 = 0.402737184e-03 * si.hPa / si.K**3
+FWC_I4 = 0.565392987e-05 * si.hPa / si.K**4
+FWC_I5 = 0.521693933e-07 * si.hPa / si.K**5
+FWC_I6 = 0.307839583e-09 * si.hPa / si.K**6
+FWC_I7 = 0.105785160e-11 * si.hPa / si.K**7
+FWC_I8 = 0.161444444e-14 * si.hPa / si.K**8
+
 rho_w = 1 * si.kilograms / si.litres
 rho_i = 916.8 * si.kg / si.metres**3
 pH_w = 7

PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py

@@ -15,3 +15,17 @@ def pvs_Celsius(T):
                 const.FWC_C7 + T *
                 const.FWC_C8
         ))))))))
+
+    @staticmethod
+    def ice_Celsius(T):
+        return (
+                const.FWC_I0 + T * (
+                const.FWC_I1 + T * (
+                const.FWC_I2 + T * (
+                const.FWC_I3 + T * (
+                const.FWC_I4 + T * (
+                const.FWC_I5 + T * (
+                const.FWC_I6 + T * (
+                const.FWC_I7 + T *
+                const.FWC_I8
+        ))))))))

PySDM/physics/spectra.py

@@ -53,6 +53,23 @@ def __init__(self, norm_factor: float, m_mode: float, s_geom: float):
     def s_geom(self):
         return math.exp(self.distribution_params[0])
 
+    @property
+    def m_mode(self):
+        return self.distribution_params[2]
+
+class TopHat:
+    def __init__(self, norm_factor, endpoints):
+        self.norm_factor = norm_factor
+        self.endpoints = endpoints
+        self._mn = endpoints[0]
+        self._mx = endpoints[1]
+
+    def cumulative(self, x):
+        return self.norm_factor * np.minimum(1, np.maximum(0, (x - self._mn) / (self._mx - self._mn)))
+
+    def percentiles(self, cdf_values):
+        return (self._mx - self._mn) * (np.asarray(cdf_values) + self._mn / (self._mx - self._mn))
+
 class Sum:
 
     def __init__(self, spectra: tuple, interpolation_grid=default_interpolation_grid):

PySDM_tests/smoke_tests/Alpert_and_Knopf_2016/test_AK16_fig_1.py

@@ -0,0 +1,87 @@
+from PySDM_examples.Alpert_and_Knopf_2016 import simulation, Table1
+from PySDM.physics import si, constants as const, Formulae
+from PySDM.physics.spectra import Lognormal
+import numpy as np
+from matplotlib import pylab
+import pytest
+
+n_runs_per_case = 3
+
+@pytest.mark.parametrize("multiplicity", (1, 2, 10))
+def test_AK16_fig_1(multiplicity, plot=False):
+    # Arrange
+    J_het = 1e3 / si.cm ** 2 / si.s
+    A_g = 1e-5 * si.cm ** 2
+
+    dt = 1 * si.s
+    total_time = 6 * si.min
+
+    # dummy multipliers (multiplied and then divided by)
+    dv = 1 * si.cm ** 3  # will become used if coalescence or other processes are turned on
+    droplet_volume = 1 * si.um ** 3  # ditto
+
+    cases = Table1(dv=dv)
+
+    # Act
+    output = {}
+
+    for key in ('Iso3', 'Iso4', 'Iso1', 'Iso2'):
+        case = cases[key]
+        output[key] = []
+        for i in range(n_runs_per_case):
+            seed = i
+            number_of_real_droplets = case['ISA'].norm_factor * dv
+            n_sd = number_of_real_droplets / multiplicity
+            assert int(n_sd) == n_sd
+            n_sd = int(n_sd)
+
+            data = simulation(seed=i, n_sd=n_sd, dt=dt, dv=dv, spectrum=case['ISA'],
+                          droplet_volume=droplet_volume, multiplicity=multiplicity, J_het=J_het,
+                          total_time=total_time, number_of_real_droplets=number_of_real_droplets)
+            output[key].append(data)
+
+    # Plot
+    if plot:
+        for key in output.keys():
+            for run in range(n_runs_per_case):
+                label = f"{key}: σ=ln({int(cases[key]['ISA'].s_geom)}),N={int(cases[key]['ISA'].norm_factor * dv)}"
+                pylab.step(
+                    dt / si.min * np.arange(len(output[key][run])),
+                    output[key][run],
+                    label=label if run == 0 else None,
+                    color=cases[key]['color'],
+                    linewidth=.666
+                )
+            output[key].append(np.mean(np.asarray(output[key]), axis=0))
+            pylab.step(
+                dt / si.min * np.arange(len(output[key][-1])),
+                output[key][-1],
+                color=cases[key]['color'],
+                linewidth=1.666
+            )
+
+        pylab.legend()
+        pylab.yscale('log')
+        pylab.ylim(1e-2, 1)
+        pylab.xlim(0, total_time / si.min)
+        pylab.xlabel("t / min")
+        pylab.ylabel("$f_{ufz}$")
+        pylab.gca().set_box_aspect(1)
+        pylab.show()
+
+    # Assert
+    np.testing.assert_array_less(
+        output['Iso3'][-1][1:int(1 * si.min / dt)],
+        output['Iso1'][-1][1:int(1 * si.min / dt)]
+    )
+    np.testing.assert_array_less(
+        output['Iso1'][-1][int(2 * si.min / dt):],
+        output['Iso3'][-1][int(2 * si.min / dt):]
+    )
+    np.testing.assert_array_less(
+        output['Iso2'][int(.5 * si.min / dt):],
+        output['Iso1'][int(.5 * si.min / dt):]
+    )
+    for key in output.keys():
+        np.testing.assert_array_less(1e-3, output[key][-1][:int(.25 * si.min / dt)])
+        np.testing.assert_array_less(output[key][-1][:], 1 + 1e-10)

PySDM_tests/unit_tests/physics/test_saturation_vapour_pressure.py

@@ -0,0 +1,44 @@
+import numpy as np
+from matplotlib import pylab
+from PySDM.physics import Formulae, constants as const
+import inspect
+
+def test_saturation_vapour_pressures(plot=False):
+    # Arrange
+    formulae = {
+        k: Formulae(saturation_vapour_pressure=k)
+        for k in ('FlatauWalkoCotton', 'AugustRocheMagnus')
+    }
+    T = np.linspace(-.2, .4)
+
+    # Plot
+    if plot:
+        pylab.axhline(const.p_tri, label='triple point', color='red')
+        pylab.axvline(const.T_tri - const.T0, color='red')
+        for k, v in formulae.items():
+            for name, func in inspect.getmembers(v.saturation_vapour_pressure):
+                if not name.startswith('__'):
+                    pylab.plot(T, func(T), label=f"{k}::{name}")
+        pylab.grid()
+        pylab.legend()
+        pylab.xlabel('T [C]')
+        pylab.ylabel('p [Pa]')
+        pylab.show()
+
+    # Assert
+    T = np.linspace(-20, 20, 100)
+    np.testing.assert_allclose(
+        Formulae(saturation_vapour_pressure='FlatauWalkoCotton').saturation_vapour_pressure.pvs_Celsius(T),
+        Formulae(saturation_vapour_pressure='AugustRocheMagnus').saturation_vapour_pressure.pvs_Celsius(T),
+        rtol=1e-2
+    )
+    T = np.linspace(-20, 0.3, 100)
+    np.testing.assert_array_less(
+        Formulae(saturation_vapour_pressure='FlatauWalkoCotton').saturation_vapour_pressure.ice_Celsius(T),
+        Formulae(saturation_vapour_pressure='FlatauWalkoCotton').saturation_vapour_pressure.pvs_Celsius(T)
+    )
+    T = np.linspace(0.35, 20, 100)
+    np.testing.assert_array_less(
+        Formulae(saturation_vapour_pressure='FlatauWalkoCotton').saturation_vapour_pressure.pvs_Celsius(T),
+        Formulae(saturation_vapour_pressure='FlatauWalkoCotton').saturation_vapour_pressure.ice_Celsius(T)
+    )

PySDM_tests/unit_tests/physics/test_spectra_top_hat.py

@@ -0,0 +1,63 @@
+from PySDM.physics import spectra
+from matplotlib import pylab
+import numpy as np
+
+
+class TestSpectraTopHat:
+    @staticmethod
+    def test_cumulative(plot=False):
+        # arrange
+        endpoints = (44, 666)
+        norm_factor = np.pi
+        spectrum = spectra.TopHat(norm_factor=norm_factor, endpoints=endpoints)
+        x = np.linspace(0, 1000)
+
+        # act
+        y = spectrum.cumulative(x)
+
+        # plot
+        if plot:
+            pylab.axhline(0)
+            pylab.axhline(norm_factor)
+            pylab.axvline(endpoints[0])
+            pylab.axvline(endpoints[1])
+            pylab.plot(x, y, color='red')
+            pylab.xlim(x[0], x[-1])
+            pylab.grid()
+            pylab.show()
+
+        # assert
+        for point in y:
+            assert 0 <= point <= norm_factor
+        assert y[-1] == norm_factor
+
+        hy, hx = np.histogram(np.diff(y) / np.diff(x))
+        assert hx[0] == 0
+        np.testing.assert_approx_equal(hx[-1], norm_factor / (endpoints[1] - endpoints[0]))
+        assert np.sum(hy[1:-1]) == 2
+
+    @staticmethod
+    def test_percentiles(plot=False):
+        # arrange
+        spectrum = spectra.TopHat(norm_factor=np.e, endpoints=(-np.pi, np.pi))
+        x = np.linspace(0, 1)
+
+        # act
+        y = spectrum.percentiles(x)
+
+        # plot
+        if plot:
+            pylab.axvline(0)
+            pylab.axvline(1)
+            pylab.axhline(spectrum.endpoints[0])
+            pylab.axhline(spectrum.endpoints[1])
+            pylab.plot(x, y, color='red')
+            pylab.show()
+
+        # assert
+        assert y[0] == spectrum.endpoints[0]
+        assert y[-1] == spectrum.endpoints[1]
+        np.testing.assert_allclose(
+            np.diff(y) / np.diff(x),
+            spectrum.endpoints[1] - spectrum.endpoints[0]
+        )

test-time-requirements.txt

@@ -4,4 +4,4 @@ ghapi
 pytest
 
 # note: if cloning both PySDM and PySDM examples, consider "pip install -e"
-PySDM-examples @ git+git://github.com/atmos-cloud-sim-uj/PySDM-examples@af1663e#egg=PySDM-examples
+PySDM-examples @ git+git://github.com/atmos-cloud-sim-uj/PySDM-examples@49a271a#egg=PySDM-examples