adding Bigg freezing spectrum; cleanups; spectral discretisation for …

…spectra without percentiles defined

PySDM/dynamics/freezing.py

@@ -1,12 +1,13 @@
-from PySDM.physics import constants
-from PySDM.physics.heterogeneous_ice_nucleation_rate import Constant, ABIFM
+from PySDM.physics.heterogeneous_ice_nucleation_rate import Null
+
 
 class Freezing:
     def __init__(self, *, singular=True):
         self.singular = singular
         self.enable = True
         self.rand = None
         self.rng = None
+        self.particulator = None
 
     def register(self, builder):
         self.particulator = builder.particulator
@@ -15,6 +16,7 @@ def register(self, builder):
         if self.singular:
             builder.request_attribute("freezing temperature")
         else:
+            assert not isinstance(builder.formulae.heterogeneous_ice_nucleation_rate, Null)
             builder.request_attribute("immersed surface area")
             self.rand = self.particulator.Storage.empty(self.particulator.n_sd, dtype=float)
             self.rng = self.particulator.Random(self.particulator.n_sd, self.particulator.bck.formulae.seed)

PySDM/initialisation/spectral_sampling.py

@@ -1,5 +1,6 @@
 import numpy as np
 from typing import Tuple
+from scipy import optimize
 from PySDM.physics import constants as const
 
 default_cdf_range = (.00001, .99999)
@@ -10,7 +11,14 @@ def __init__(self, spectrum, size_range: [None, Tuple[float, float]] = None):
         self.spectrum = spectrum
 
         if size_range is None:
-            self.size_range = spectrum.percentiles(default_cdf_range)
+            if hasattr(spectrum, 'percentiles'):
+                self.size_range = spectrum.percentiles(default_cdf_range)
+            else:
+                self.size_range = [np.nan, np.nan]
+                for i in (0, 1):
+                    result = optimize.root(lambda x: spectrum.cdf(x) - default_cdf_range[i], x0=spectrum.median())
+                    assert result.success
+                    self.size_range[i] = result.x
         else:
             assert len(size_range) == 2
             assert size_range[0] > 0
@@ -69,12 +77,13 @@ def sample(self, n_sd):
 
         return self._sample(percentiles, self.spectrum)
 
+
 class UniformRandom(SpectralSampling):
     def __init__(self, spectrum, size_range=None, seed=const.default_random_seed):
         super().__init__(spectrum, size_range)
         self.rng = np.random.default_rng(seed)
 
     def sample(self, n_sd):
         pdf_arg = self.rng.uniform(*self.size_range, n_sd)
-        dr = (self.size_range[1] - self.size_range[0]) / n_sd
+        dr = abs(self.size_range[1] - self.size_range[0]) / n_sd
         return pdf_arg, dr * self.spectrum.size_distribution(pdf_arg)

PySDM/physics/freezing_temperature_spectrum/__init__.py

@@ -1,2 +1,3 @@
 from .niemand_et_al_2012 import Niemand_et_al_2012
+from .bigg_1953 import Bigg_1953
 from .null import Null

PySDM/physics/freezing_temperature_spectrum/bigg_1953.py

@@ -0,0 +1,27 @@
+import numpy as np
+from PySDM.physics import constants as const
+
+P_median = .5
+DT_median = np.nan
+
+# TODO #599: there are two Bigg 1953 papers
+# TODO #599: relate DT to drop volume to A_insol? (the second paper!)
+
+
+class Bigg_1953:
+    def __init__(self):
+        assert np.isfinite(DT_median)
+
+    @staticmethod
+    def pdf(T, A_insol):
+        A = np.log(1 - P_median)
+        B = DT_median - const.T0
+        return - A * np.exp(A * np.exp(B + T) + B + T)
+
+    @staticmethod
+    def cdf(T, A_insol):
+        return np.exp(np.log(1 - P_median) * np.exp(DT_median - (const.T0 - T)))
+
+    @staticmethod
+    def median():
+        return const.T0 - DT_median

PySDM/physics/freezing_temperature_spectrum/niemand_et_al_2012.py

@@ -1,20 +1,24 @@
-from PySDM.physics import constants as const, si
+from PySDM.physics import constants as const
 import numpy as np
 
-a = -0.517
-b = 8.934
-A_insol = const.pi * (1*si.um)**2
+a = np.nan
+b = np.nan
+
 
 class Niemand_et_al_2012:
     def __str__(self):
         return 'Niemand et al. 2012'
 
+    def __init__(self):
+        assert np.isfinite(a)
+        assert np.isfinite(b)
+
     @staticmethod
-    def pdf(T):
+    def pdf(T, A_insol):
         ns_T = np.exp(a * (T - const.T0) + b)
         return -A_insol * a * ns_T * np.exp(-A_insol * ns_T)
 
     @staticmethod
-    def cdf(T):
+    def cdf(T, A_insol):
         ns_T = np.exp(a * (T - const.T0) + b)
         return 1 - np.exp(-A_insol * ns_T) - np.exp(-A_insol*np.exp(-a * const.T0 + b))

PySDM/physics/freezing_temperature_spectrum/null.py

@@ -1,2 +1,4 @@
 class Null:
-    pass
+    @staticmethod
+    def cdf(T):
+        pass

PySDM/physics/heterogeneous_ice_nucleation_rate/abifm.py

@@ -5,6 +5,7 @@
 c = np.inf
 unit = 1 / si.cm**2 / si.s
 
+
 class ABIFM:
     @staticmethod
     def _check():

PySDM/physics/spectra.py

@@ -22,6 +22,9 @@ def size_distribution(self, x):
         result = self.norm_factor * self.distribution.pdf(x, *self.distribution_params)
         return result
 
+    def pdf(self, x):
+        return self.size_distribution(x) / self.norm_factor
+
     def stats(self, moments):
         result = self.distribution.stats(*self.distribution_params, moments)
         return result
@@ -57,6 +60,10 @@ def s_geom(self):
     def m_mode(self):
         return self.distribution_params[2]
 
+    def __str__(self):
+        return f"{self.__class__.__name__}:(N={self.norm_factor}, m_mode={self.m_mode}, s_geom={self.s_geom})"
+
+
 class TopHat:
     def __init__(self, norm_factor, endpoints):
         self.norm_factor = norm_factor
@@ -70,6 +77,7 @@ def cumulative(self, x):
     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):

tests/unit_tests/initialisation/test_spectral_discretisation.py

@@ -1,33 +1,35 @@
 from PySDM.initialisation import spectral_sampling, spectro_glacial
 from PySDM.physics.spectra import Lognormal
 from PySDM.physics import Formulae, constants as const
+from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012
 import numpy as np
 import pytest
 
+niemand_et_al_2012.a = -0.517
+niemand_et_al_2012.b = 8.934
+
 m_mode = .5e-5
 n_part = 256 * 16
 s_geom = 1.5
 spectrum = Lognormal(n_part, m_mode, s_geom)
 m_range = (.1e-6, 100e-6)
-formulae = Formulae()
+formulae = Formulae(freezing_temperature_spectrum='Niemand_et_al_2012', seed=const.default_random_seed)
+
 
 @pytest.mark.parametrize("discretisation", [
 	pytest.param(spectral_sampling.Linear(spectrum, m_range)),
 	pytest.param(spectral_sampling.Logarithmic(spectrum, m_range)),
 	pytest.param(spectral_sampling.ConstantMultiplicity(spectrum, m_range)),
 	pytest.param(spectral_sampling.UniformRandom(spectrum, m_range)),
-	pytest.param(spectro_glacial.Independent(
-		size_spectrum=spectrum,
-		freezing_temperature_spectrum=formulae.freezing_temperature_spectrum
-	))
+	# TODO #599
 ])
 def test_spectral_discretisation(discretisation):
 	# Arrange
-	n_sd = 10000
+	n_sd = 100000
 
 	# Act
 	if isinstance(discretisation, spectro_glacial.SpectroGlacialSampling):
-		m, _, n = discretisation.sample(n_sd)
+		m, _, __, n = discretisation.sample(n_sd)
 	else:
 		m, n = discretisation.sample(n_sd)
 

tests/unit_tests/physics/test_freezing_temperature_spectra.py

@@ -1,32 +1,47 @@
-from PySDM.physics import Formulae, constants as const
-from matplotlib import pylab
 import numpy as np
+import pytest
+from matplotlib import pylab
+from PySDM.physics import Formulae, si, constants as const
+from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012, bigg_1953
+
 
+A = 1 * si.um**2
 
-def test_freezing_temperature_spectra(plot=False):
+
+@pytest.mark.parametrize("model", (
+        'Niemand_et_al_2012',
+        'Bigg_1953'
+))
+def test_freezing_temperature_spectra(model, plot=False):
     # Arrange
-    formulae = Formulae()
+    bigg_1953.DT_median = 33
+    niemand_et_al_2012.a = -0.517
+    niemand_et_al_2012.b = 8.934
+    niemand_et_al_2012.A_insol = 1 * si.um ** 2
+
+    formulae = Formulae(freezing_temperature_spectrum=model)
     T = np.linspace(const.T0, const.T0 - 40, num=100)
 
     # Act
-    pdf = formulae.freezing_temperature_spectrum.pdf(T)
-    cdf = formulae.freezing_temperature_spectrum.cdf(T)
+    pdf = formulae.freezing_temperature_spectrum.pdf(T, A)
+    cdf = formulae.freezing_temperature_spectrum.cdf(T, A)
 
     # Plot
+    pylab.plot(T, pdf, linestyle='-', marker='o', label='pdf')
+    pdfax = pylab.gca()
+    cdfax = pdfax.twinx()
+    cdfax.plot(T, cdf, linestyle='--', marker='x', label='cdf')
+    pylab.xlabel('T [K]')
+    pylab.xlim(np.amax(T), np.amin(T))
+    pdfax.set_ylabel('pdf [K$^{-1}$]')
+    cdfax.set_ylabel('cdf [1]')
+    pylab.grid()
+    pylab.title(model)
     if plot:
-        pylab.plot(T, pdf, linestyle='-', marker='o', label='pdf')
-        pdfax = pylab.gca()
-        cdfax = pdfax.twinx()
-        cdfax.plot(T, cdf, linestyle='--', marker='x', label='cdf')
-        pylab.xlabel('T [K]')
-        pylab.xlim(np.amax(T), np.amin(T))
-        pdfax.set_ylabel('pdf [K$^{-1}$]')
-        cdfax.set_ylabel('cdf [1]')
-        pylab.grid()
         pylab.show()
 
     # Assert
     dT = abs(T[1] - T[0])
-    np.testing.assert_approx_equal(np.sum(pdf * dT), 1)
-    np.testing.assert_approx_equal(cdf[0]+1, 1)
-    np.testing.assert_approx_equal(cdf[-1], 1)
+    np.testing.assert_approx_equal(np.sum(pdf * dT), 1, significant=3)
+    np.testing.assert_approx_equal(cdf[0]+1, 1, significant=3)
+    np.testing.assert_approx_equal(cdf[-1], 1, significant=3)