Merge pull request #658 from slayoo/freezing

cleanups, null default for freezing spectrum, README update

Author: slayoo

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():