Merge pull request #658 from slayoo/freezing

cleanups, null default for freezing spectrum, README update

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/environments/box.py

@@ -18,7 +18,10 @@ def __getitem__(self, item):
         return self._ambient_air[item]
 
     def __setitem__(self, key, value):
-        self._ambient_air[key] = self.particulator.backend.Storage.from_ndarray(np.array([value]))
+        if key not in self._ambient_air:
+            self._ambient_air[key] = self.particulator.backend.Storage.from_ndarray(np.array([value]))
+        else:
+            self._ambient_air[key][:] = value
 
     def register(self, builder):
         self.particulator = builder.particulator

PySDM/initialisation/multiplicities.py

@@ -12,6 +12,6 @@ def discretise_n(y_float):
         raise Exception(f"{percent_diff}% error in total real-droplet number due to casting multiplicities to ints")
 
     if not (y_int > 0).all():
-        raise Exception("int-casting resulted in multiplicity of zero")
+        raise Exception(f"int-casting resulted in multiplicity of zero (min(y_float)={min(y_float)})")
 
     return y_int

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/formulae.py

@@ -65,7 +65,7 @@ def _pick(value: str, choices: dict):
     for name, cls in choices.items():
         if name == value:
             return cls()
-    raise ValueError(f"Unknown setting: '{value}';, choices are: {tuple(choices.keys())}")
+    raise ValueError(f"Unknown setting: '{value}'; choices are: {tuple(choices.keys())}")
 
 
 def _choices(module):
@@ -93,7 +93,7 @@ def __init__(self, *,
                  state_variable_triplet: str = 'RhodThdQv',
                  particle_advection: str = 'ImplicitInSpace',
                  hydrostatics: str = 'Default',
-                 freezing_temperature_spectrum: str = 'Niemand_et_al_2012',
+                 freezing_temperature_spectrum: str = 'Null',
                  heterogeneous_ice_nucleation_rate: str = 'Null'
                  ):
         self.seed = seed

PySDM/physics/freezing_temperature_spectrum/__init__.py

@@ -1 +1,3 @@
-from .niemand_et_al_2012 import Niemand_et_al_2012
+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

@@ -0,0 +1,4 @@
+class Null:
+    @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):

PySDM/products/product.py

@@ -79,14 +79,14 @@ def __init__(self, **args):
         self.source = None
 
     def register(self, builder):
-        assert isinstance(builder.particulator.env, _Moist)
         super().register(builder)
         self.particulator.observers.append(self)
         self.environment = builder.particulator.env
         self.source = self.environment[self._name]
 
     def notify(self):
-        self.source = self.environment.get_predicted(self._name)
+        if isinstance(self.environment, _Moist):
+            self.source = self.environment.get_predicted(self._name)
 
     def get(self):
         self.download_to_buffer(self.source)

PySDM/state/particles.py

@@ -76,6 +76,9 @@ def keys(self):
     def __getitem__(self, item):
         return self.attributes[item].get()
 
+    def __contains__(self, key):
+        return key in self.attributes
+
     def permutation(self, u01, local):
         if local:
             """

README.md

@@ -602,7 +602,8 @@ the license of the contributed code must be compatible with GPL v3.
 
 Developing the code, we follow [The Way of Python](https://www.python.org/dev/peps/pep-0020/) and 
 the [KISS principle](https://en.wikipedia.org/wiki/KISS_principle).
-The codebase has greatly benefited from [PyCharm code inspections](https://www.jetbrains.com/help/pycharm/code-inspection.html).
+The codebase has greatly benefited from [PyCharm code inspections](https://www.jetbrains.com/help/pycharm/code-inspection.html)
+and [Pylint](https://pylint.org) code analysis (which constitutes one of the CI workflows).
 
 Issues regarding any incorrect, unintuitive or undocumented bahaviour of
 PySDM are best to be reported on the [GitHub issue tracker](https://github.com/atmos-cloud-sim-uj/PySDM/issues/new).
@@ -616,10 +617,24 @@ We look forward to your contributions and feedback.
 
 ## Credits:
 
-Development of PySDM is supported by the EU through a grant of the Foundation for Polish Science (POIR.04.04.00-00-5E1C/18).
-
-copyright: Jagiellonian University   
-licence: GPL v3   
+The development and maintenance of PySDM is led by [Sylwester Arabas](https://github.com/slayoo/).
+[Piotr Bartman](https://github.com/piotrbartman/) had been the architect and main developer 
+of technological solutions in PySDM. 
+The suite of examples shipped with PySDM includes contributions from researchers 
+from [Jagiellonian University](https://en.uj.edu.pl/en) departments of computer science, physics and chemistry;
+and from 
+[Caltech's Climate Modelling Alliance](https://clima.caltech.edu/).
+
+Development of PySDM had been initially supported by the EU through a grant of the 
+[Foundation for Polish Science](https://www.fnp.org.pl/)) (POIR.04.04.00-00-5E1C/18) 
+realised at the [Jagiellonian University](https://en.uj.edu.pl/en).
+The immersion freezing support in PySDM is developed with support from the
+US Department of Energy [Atmospheric System Research](https://asr.science.energy.gov/) programme
+through a grant realised at the 
+[University of Illinois at Urbana-Champaign](https://illinois.edu/).
+
+copyright: [Jagiellonian University](https://en.uj.edu.pl/en)    
+licence: [GPL v3](https://www.gnu.org/licenses/gpl-3.0.html)
 
 ## Related resources and open-source projects
 

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)