new tests for mixed-phase moments, signed-water-mass attribute, radius_to_mass() for mixed phased spheres (#1513)

Co-authored-by: Sylwester Arabas <sylwester.arabas@agh.edu.pl>

Author: tluettm

PySDM/attributes/isotopes/moles.py

@@ -62,7 +62,7 @@ def __init__(self, builder):
                 / M_H2O,
                 builder.get_attribute("moles_18O"): -(2 * const.M_1H + const.M_18O)
                 / M_H2O,
-                builder.get_attribute("water mass"): 1 / M_H2O,
+                builder.get_attribute("signed water mass"): 1 / M_H2O,
             },
         )
 

PySDM/attributes/physics/__init__.py

@@ -10,7 +10,7 @@
 from .equilibrium_supersaturation import EquilibriumSupersaturation
 from .heat import Heat
 from .hygroscopicity import Kappa, KappaTimesDryVolume
-from .water_mass import WaterMass
+from .water_mass import SignedWaterMass
 from .multiplicity import Multiplicity
 from .radius import Radius, SquareRootOfRadius
 from .relative_fall_velocity import RelativeFallMomentum, RelativeFallVelocity

PySDM/attributes/physics/water_mass.py

@@ -1,13 +1,63 @@
 """
-particle mass, derived attribute for coalescence
-in simulation involving mixed-phase clouds, positive values correspond to
+particle mass attributes
+in simulations involving mixed-phase clouds, positive values correspond to
 liquid water and negative values to ice
 """
 
-from PySDM.attributes.impl import ExtensiveAttribute, register_attribute
+from PySDM.attributes.impl import (
+    ExtensiveAttribute,
+    DerivedAttribute,
+    register_attribute,
+)
 
 
 @register_attribute()
-class WaterMass(ExtensiveAttribute):
+class SignedWaterMass(ExtensiveAttribute):
     def __init__(self, builder):
-        super().__init__(builder, name="water mass")
+        super().__init__(builder, name="signed water mass")
+
+
+@register_attribute(
+    name="water mass",
+    variant=lambda _, formulae: not formulae.particle_shape_and_density.supports_mixed_phase(),
+)
+class ViewWaterMass(DerivedAttribute):
+    def __init__(self, builder):
+        self.signed_water_mass = builder.get_attribute("signed water mass")
+
+        super().__init__(
+            builder,
+            name="water mass",
+            dependencies=(self.signed_water_mass,),
+        )
+
+    def mark_updated(self):
+        self.signed_water_mass.mark_updated()
+
+    def allocate(self, idx):
+        pass
+
+    def recalculate(self):
+        pass
+
+    def get(self):
+        return self.signed_water_mass.data
+
+
+@register_attribute(
+    name="water mass",
+    variant=lambda _, formulae: formulae.particle_shape_and_density.supports_mixed_phase(),
+)
+class AbsWaterMass(DerivedAttribute):
+    def __init__(self, builder):
+        self.signed_water_mass = builder.get_attribute("signed water mass")
+
+        super().__init__(
+            builder,
+            name="water mass",
+            dependencies=(self.signed_water_mass,),
+        )
+
+    def recalculate(self):
+        self.data.fill(self.signed_water_mass.data)
+        self.data.abs()

PySDM/backends/impl_common/freezing_attributes.py

@@ -1,20 +1,26 @@
 """
-groups of attributes used in either singular or time-dependent immmersion freezing regimes
+groups of attributes used in either singular or time-dependent immersion freezing regimes
 """
 
 from collections import namedtuple
 
 
 class SingularAttributes(
-    namedtuple("SingularAttributes", ("freezing_temperature", "water_mass"))
+    namedtuple(
+        typename="SingularAttributes",
+        field_names=("freezing_temperature", "signed_water_mass"),
+    )
 ):
     """groups attributes required in singular regime"""
 
     __slots__ = ()
 
 
 class TimeDependentAttributes(
-    namedtuple("TimeDependentAttributes", ("immersed_surface_area", "water_mass"))
+    namedtuple(
+        typename="TimeDependentAttributes",
+        field_names=("immersed_surface_area", "signed_water_mass"),
+    )
 ):
     """groups attributes required in time-dependent regime"""
 

PySDM/backends/impl_numba/methods/freezing_methods.py

@@ -39,16 +39,16 @@ def freeze_singular_body(
                 if attributes.freezing_temperature[i] == 0:
                     continue
                 if thaw and frozen_and_above_freezing_point(
-                    attributes.water_mass[i], temperature[cell[i]]
+                    attributes.signed_water_mass[i], temperature[cell[i]]
                 ):
-                    _thaw(attributes.water_mass, i)
+                    _thaw(attributes.signed_water_mass, i)
                 elif (
                     unfrozen_and_saturated(
-                        attributes.water_mass[i], relative_humidity[cell[i]]
+                        attributes.signed_water_mass[i], relative_humidity[cell[i]]
                     )
                     and temperature[cell[i]] <= attributes.freezing_temperature[i]
                 ):
-                    _freeze(attributes.water_mass, i)
+                    _freeze(attributes.signed_water_mass, i)
 
         self.freeze_singular_body = freeze_singular_body
 
@@ -68,17 +68,17 @@ def freeze_time_dependent_body(  # pylint: disable=unused-argument,too-many-argu
             freezing_temperature,
             thaw,
         ):
-            n_sd = len(attributes.water_mass)
+            n_sd = len(attributes.signed_water_mass)
             for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
                 if attributes.immersed_surface_area[i] == 0:
                     continue
                 cell_id = cell[i]
                 if thaw and frozen_and_above_freezing_point(
-                    attributes.water_mass[i], temperature[cell_id]
+                    attributes.signed_water_mass[i], temperature[cell_id]
                 ):
-                    _thaw(attributes.water_mass, i)
+                    _thaw(attributes.signed_water_mass, i)
                 elif unfrozen_and_saturated(
-                    attributes.water_mass[i], relative_humidity[cell_id]
+                    attributes.signed_water_mass[i], relative_humidity[cell_id]
                 ):
                     rate_assuming_constant_temperature_within_dt = (
                         j_het(a_w_ice[cell_id]) * attributes.immersed_surface_area[i]
@@ -87,7 +87,7 @@ def freeze_time_dependent_body(  # pylint: disable=unused-argument,too-many-argu
                         r=rate_assuming_constant_temperature_within_dt, dt=timestep
                     )
                     if rand[i] < prob:
-                        _freeze(attributes.water_mass, i)
+                        _freeze(attributes.signed_water_mass, i)
                         # if record_freezing_temperature:
                         #     freezing_temperature[i] = temperature[cell_id]
 
@@ -99,7 +99,7 @@ def freeze_singular(
         self.freeze_singular_body(
             SingularAttributes(
                 freezing_temperature=attributes.freezing_temperature.data,
-                water_mass=attributes.water_mass.data,
+                signed_water_mass=attributes.signed_water_mass.data,
             ),
             temperature.data,
             relative_humidity.data,
@@ -125,7 +125,7 @@ def freeze_time_dependent(
             rand.data,
             TimeDependentAttributes(
                 immersed_surface_area=attributes.immersed_surface_area.data,
-                water_mass=attributes.water_mass.data,
+                signed_water_mass=attributes.signed_water_mass.data,
             ),
             timestep,
             cell.data,

PySDM/backends/impl_numba/storage.py

@@ -208,3 +208,6 @@ def fill(self, other):
 
     def exp(self):
         self.data[:] = np.exp(self.data)
+
+    def abs(self):
+        self.data[:] = np.abs(self.data)

PySDM/backends/impl_thrust_rtc/methods/freezing_methods.py

@@ -18,7 +18,7 @@ def freeze_time_dependent_body(self):
             param_names=(
                 "rand",
                 "immersed_surface_area",
-                "water_mass",
+                "signed_water_mass",
                 "timestep",
                 "cell",
                 "a_w_ice",
@@ -32,12 +32,12 @@ def freeze_time_dependent_body(self):
                     return;
                 }}
                 if (thaw && {self.formulae.trivia.frozen_and_above_freezing_point.c_inline(
-                    water_mass="water_mass[i]",
+                    signed_water_mass="signed_water_mass[i]",
                     temperature="temperature[cell[i]]"
                 )}) {{
-                    water_mass[i] = -1 * water_mass[i];
+                    signed_water_mass[i] = -1 * signed_water_mass[i];
                 }} else if ({self.formulae.trivia.unfrozen_and_saturated.c_inline(
-                        water_mass="water_mass[i]",
+                        signed_water_mass="signed_water_mass[i]",
                         relative_humidity="relative_humidity[cell[i]]"
                     )}) {{
                     auto rate_assuming_constant_temperature_within_dt = {self.formulae.heterogeneous_ice_nucleation_rate.j_het.c_inline(
@@ -48,7 +48,7 @@ def freeze_time_dependent_body(self):
                         dt="timestep"
                     )};
                     if (rand[i] < prob) {{
-                        water_mass[i] = -1 * water_mass[i];
+                        signed_water_mass[i] = -1 * signed_water_mass[i];
                     }}
                 }}
             """.replace(
@@ -61,7 +61,7 @@ def freeze_singular_body(self):
         return trtc.For(
             param_names=(
                 "freezing_temperature",
-                "water_mass",
+                "signed_water_mass",
                 "temperature",
                 "relative_humidity",
                 "cell",
@@ -73,17 +73,17 @@ def freeze_singular_body(self):
                     return;
                 }}
                 if (thaw && {self.formulae.trivia.frozen_and_above_freezing_point.c_inline(
-                    water_mass="water_mass[i]",
+                    signed_water_mass="signed_water_mass[i]",
                     temperature="temperature[cell[i]]"
                 )}) {{
-                    water_mass[i] = -1 * water_mass[i];
+                    signed_water_mass[i] = -1 * signed_water_mass[i];
                 }} else if (
                     {self.formulae.trivia.unfrozen_and_saturated.c_inline(
-                        water_mass="water_mass[i]",
+                        signed_water_mass="signed_water_mass[i]",
                         relative_humidity="relative_humidity[cell[i]]"
                     )} && temperature[cell[i]] <= freezing_temperature[i]
                 ) {{
-                    water_mass[i] = -1 * water_mass[i];
+                    signed_water_mass[i] = -1 * signed_water_mass[i];
                 }}
             """.replace(
                 "real_type", self._get_c_type()
@@ -99,7 +99,7 @@ def freeze_singular(
             n=n_sd,
             args=(
                 attributes.freezing_temperature.data,
-                attributes.water_mass.data,
+                attributes.signed_water_mass.data,
                 temperature.data,
                 relative_humidity.data,
                 cell.data,
@@ -128,7 +128,7 @@ def freeze_time_dependent(  # pylint: disable=unused-argument
             args=(
                 rand.data,
                 attributes.immersed_surface_area.data,
-                attributes.water_mass.data,
+                attributes.signed_water_mass.data,
                 self._get_floating_point(timestep),
                 cell.data,
                 a_w_ice.data,

PySDM/backends/impl_thrust_rtc/storage.py

@@ -319,6 +319,19 @@ def divide_if_not_zero(output, divisor):
         def exp(output):
             Impl.__exp_body.launch_n(output.shape[0], Impl.thrust((output,)))
 
+        __abs_body = trtc.For(
+            ("output",),
+            "i",
+            """
+                output[i] = abs(output[i]);
+            """,
+        )
+
+        @staticmethod
+        @nice_thrust(**NICE_THRUST_FLAGS)
+        def abs(output):
+            Impl.__abs_body.launch_n(output.shape[0], Impl.thrust((output,)))
+
     class Storage(StorageBase):
         FLOAT = BACKEND._get_np_dtype()
         INT = np.int64
@@ -554,4 +567,8 @@ def exp(self):
             Impl.exp(self)
             return self
 
+        def abs(self):
+            Impl.abs(self)
+            return self
+
     return Storage

PySDM/builder.py

@@ -109,12 +109,19 @@ def build(
             ), "implied volume-to-mass conversion is only supported for spherical particles"
             attributes["water mass"] = (
                 self.particulator.formulae.particle_shape_and_density.volume_to_mass(
-                    attributes["volume"]
+                    attributes.pop("volume")
                 )
             )
-            del attributes["volume"]
             self.request_attribute("volume")
 
+        if (
+            "water mass" in attributes
+            and "signed water mass" not in attributes
+            and not self.particulator.formulae.particle_shape_and_density.supports_mixed_phase()
+        ):
+            attributes["signed water mass"] = attributes.pop("water mass")
+            self.request_attribute("water mass")
+
         self.req_attr = {}
         for attr_name in self.req_attr_names:
             self._resolve_attribute(attr_name)

PySDM/dynamics/freezing.py

@@ -27,7 +27,7 @@ def register(self, builder):
 
         assert builder.formulae.particle_shape_and_density.supports_mixed_phase()
 
-        builder.request_attribute("volume")
+        builder.request_attribute("signed water mass")
         if self.singular or self.record_freezing_temperature:
             builder.request_attribute("freezing temperature")
 
@@ -59,7 +59,7 @@ def __call__(self):
                     freezing_temperature=self.particulator.attributes[
                         "freezing temperature"
                     ],
-                    water_mass=self.particulator.attributes["water mass"],
+                    signed_water_mass=self.particulator.attributes["signed water mass"],
                 ),
                 temperature=self.particulator.environment["T"],
                 relative_humidity=self.particulator.environment["RH"],
@@ -74,7 +74,7 @@ def __call__(self):
                     immersed_surface_area=self.particulator.attributes[
                         "immersed surface area"
                     ],
-                    water_mass=self.particulator.attributes["water mass"],
+                    signed_water_mass=self.particulator.attributes["signed water mass"],
                 ),
                 timestep=self.particulator.dt,
                 cell=self.particulator.attributes["cell id"],
@@ -90,6 +90,6 @@ def __call__(self):
                 thaw=self.thaw,
             )
 
-        self.particulator.attributes.mark_updated("water mass")
+        self.particulator.attributes.mark_updated("signed water mass")
         if self.record_freezing_temperature:
             self.particulator.attributes.mark_updated("freezing temperature")