improve test coverage for Mesh; introduce Mesh.domain_bottom_surface_area; switch from volume to mass for rainfall counting; improve docstrings (default length of omitted dims, liquid/ice phase). closes #708

PySDM/backends/impl_numba/methods/displacement_methods.py

@@ -152,15 +152,15 @@ def _flag_precipitated_body(self):
         def body(
             cell_origin,
             position_in_cell,
-            volume,
+            water_mass,
             multiplicity,
             idx,
             length,
             healthy,
             precipitation_counting_level_index,
             displacement,
         ):
-            rainfall = 0.0
+            rainfall_mass = 0.0
             flag = len(idx)
             for i in range(length):
                 position_within_column = (
@@ -173,10 +173,10 @@ def body(
                     # and crossed precip-counting level
                     position_within_column < precipitation_counting_level_index
                 ):
-                    rainfall += volume[idx[i]] * multiplicity[idx[i]]  # TODO #599
+                    rainfall_mass += abs(water_mass[idx[i]]) * multiplicity[idx[i]]
                     idx[i] = flag
                     healthy[0] = 0
-            return rainfall
+            return rainfall_mass
 
         return body
 
@@ -204,20 +204,23 @@ def body(
     # pylint: disable=too-many-arguments
     def flag_precipitated(
         self,
+        *,
         cell_origin,
         position_in_cell,
-        volume,
+        water_mass,
         multiplicity,
         idx,
         length,
         healthy,
         precipitation_counting_level_index,
         displacement,
     ) -> float:
+        """return a scalar value corresponding to the mass of water (all phases) that crossed
+        the bottom boundary of the entire domain"""
         return self._flag_precipitated_body(
             cell_origin.data,
             position_in_cell.data,
-            volume.data,
+            water_mass.data,
             multiplicity.data,
             idx.data,
             length,

PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py

@@ -73,16 +73,19 @@ def __flag_precipitated_body(self):
                 "healthy",
                 "cell_origin",
                 "position_in_cell",
-                "volume",
+                "water_mass",
                 "multiplicity",
-                "rainfall",
+                "rainfall_mass",
             ),
             "i",
             """
             auto origin = cell_origin[n_sd * (n_dims-1) + idx[i]];
             auto pic = position_in_cell[n_sd * (n_dims-1) + idx[i]];
             if (origin + pic < 0) {
-                atomicAdd((real_type*) &rainfall[0], multiplicity[idx[i]] * volume[idx[i]]);
+                atomicAdd(
+                    (real_type*) &rainfall_mass[0],
+                    multiplicity[idx[i]] * abs(water_mass[idx[i]])
+                );
                 idx[i] = n_sd;
                 healthy[0] = 0;
             }
@@ -118,7 +121,7 @@ def flag_precipitated(  # pylint: disable=unused-argument
         *,
         cell_origin,
         position_in_cell,
-        volume,
+        water_mass,
         multiplicity,
         idx,
         length,
@@ -130,8 +133,8 @@ def flag_precipitated(  # pylint: disable=unused-argument
             raise NotImplementedError()
         n_sd = trtc.DVInt64(cell_origin.shape[1])
         n_dims = trtc.DVInt64(len(cell_origin.shape))
-        rainfall = trtc.device_vector(self._get_c_type(), 1)
-        trtc.Fill(rainfall, self._get_floating_point(0))
+        rainfall_mass = trtc.device_vector(self._get_c_type(), 1)
+        trtc.Fill(rainfall_mass, self._get_floating_point(0))
         self.__flag_precipitated_body.launch_n(
             length,
             (
@@ -141,12 +144,12 @@ def flag_precipitated(  # pylint: disable=unused-argument
                 healthy.data,
                 cell_origin.data,
                 position_in_cell.data,
-                volume.data,
+                water_mass.data,
                 multiplicity.data,
-                rainfall,
+                rainfall_mass,
             ),
         )
-        return rainfall.to_host()[0]
+        return rainfall_mass.to_host()[0]
 
     @staticmethod
     def flag_out_of_column(  # pylint: disable=unused-argument

PySDM/dynamics/displacement.py

@@ -32,7 +32,7 @@ def __init__(
         self.courant = None
         self.displacement = None
         self.temp = None
-        self.precipitation_in_last_step = 0
+        self.precipitation_mass_in_last_step = 0
         self.precipitation_counting_level_index = precipitation_counting_level_index
 
         self.adaptive = adaptive
@@ -102,14 +102,14 @@ def __call__(self):
         cell_origin = self.particulator.attributes["cell origin"]
         position_in_cell = self.particulator.attributes["position in cell"]
 
-        self.precipitation_in_last_step = 0.0
+        self.precipitation_mass_in_last_step = 0.0
         for _ in range(self._n_substeps):
             self.calculate_displacement(
                 self.displacement, self.courant, cell_origin, position_in_cell
             )
             self.update_position(position_in_cell, self.displacement)
             if self.enable_sedimentation:
-                self.precipitation_in_last_step += self.particulator.remove_precipitated(
+                self.precipitation_mass_in_last_step += self.particulator.remove_precipitated(
                     displacement=self.displacement,
                     precipitation_counting_level_index=self.precipitation_counting_level_index,
                 )

PySDM/impl/mesh.py

@@ -3,11 +3,13 @@
 """
 
 import numpy as np
+from PySDM.physics import si
 
 
 class Mesh:
     # pylint: disable=too-many-arguments
     def __init__(self, grid, size, n_cell=None, dv=None, n_dims=None, strides=None):
+        """sizes of dimensions not specified in grid/size are assumed to be of 1 m"""
         self.grid = grid
         self.size = size
         self.strides = strides or Mesh.__strides(grid)
@@ -27,6 +29,14 @@ def dimension(self):
     def dim(self):
         return self.n_dims
 
+    @property
+    def domain_bottom_surface_area(self):
+        assert self.n_dims > 0
+        return {
+            1: 1 * si.m**2,
+            2: 1 * si.m * self.size[0],
+        }[self.n_dims]
+
     @staticmethod
     def mesh_0d(dv=None):
         return Mesh(

PySDM/particulator.py

@@ -396,10 +396,10 @@ def adaptive_sdm_end(self, dt_left):
     def remove_precipitated(
         self, *, displacement, precipitation_counting_level_index
     ) -> float:
-        res = self.backend.flag_precipitated(
+        rainfall_mass = self.backend.flag_precipitated(
             cell_origin=self.attributes["cell origin"],
             position_in_cell=self.attributes["position in cell"],
-            volume=self.attributes["volume"],
+            water_mass=self.attributes["water mass"],
             multiplicity=self.attributes["multiplicity"],
             idx=self.attributes._ParticleAttributes__idx,
             length=self.attributes.super_droplet_count,
@@ -408,7 +408,7 @@ def remove_precipitated(
             displacement=displacement,
         )
         self.attributes.sanitize()
-        return res
+        return rainfall_mass
 
     def flag_out_of_column(self):
         self.backend.flag_out_of_column(

PySDM/products/displacement/surface_precipitation.py

@@ -1,5 +1,6 @@
 """
-water volume flux derived from sizes of particles crossing bottom domain boundary
+water (equivalent) volume flux derived from sizes of particles crossing bottom domain boundary
+(computed from mass of both liquid and ice water)
 """
 
 from PySDM.products.impl import Product, register_product
@@ -10,36 +11,37 @@
     def __init__(self, name=None, unit="m/s"):
         super().__init__(unit=unit, name=name)
         self.displacement = None
-        self.dv = None
-        self.dz = None
+        self.domain_bottom_surface_area = None
         self._reset_counters()
 
     def _reset_counters(self):
-        self.accumulated_rainfall = 0.0
+        self.accumulated_rainfall_mass = 0.0
         self.elapsed_time = 0.0
 
     def register(self, builder):
         super().register(builder)
         self.particulator.observers.append(self)
         self.shape = ()
         self.displacement = self.particulator.dynamics["Displacement"]
-        self.dv = self.particulator.mesh.dv
-        self.dz = self.particulator.mesh.dz
+        self.domain_bottom_surface_area = (
+            self.particulator.mesh.domain_bottom_surface_area
+        )
 
     def _impl(self, **kwargs) -> float:
         if self.elapsed_time == 0.0:
             return 0.0
 
-        # TODO #708
         result = (
-            self.formulae.constants.rho_w
-            * self.accumulated_rainfall
+            self.accumulated_rainfall_mass
+            / self.formulae.constants.rho_w
             / self.elapsed_time
-            / (self.dv / self.dz)
+            / self.domain_bottom_surface_area
         )
         self._reset_counters()
         return result
 
     def notify(self):
-        self.accumulated_rainfall += self.displacement.precipitation_in_last_step
+        self.accumulated_rainfall_mass += (
+            self.displacement.precipitation_mass_in_last_step
+        )
         self.elapsed_time += self.displacement.particulator.dt

tests/unit_tests/dynamics/displacement/test_sedimentation.py

@@ -26,12 +26,12 @@ def test_boundary_condition(backend_class):
         particulator.attributes._ParticleAttributes__attributes[
             "relative fall velocity"
         ] = ConstantTerminalVelocity(particulator.backend, particulator)
-        assert sut.precipitation_in_last_step == 0
+        assert sut.precipitation_mass_in_last_step == 0
 
         # Act
         sut()
         particulator.attributes.sanitize()
 
         # Assert
         assert particulator.attributes.super_droplet_count == 0
-        assert sut.precipitation_in_last_step != 0
+        assert sut.precipitation_mass_in_last_step != 0

tests/unit_tests/impl/test_mesh.py

@@ -3,6 +3,7 @@
 import pytest
 
 from PySDM.impl.mesh import Mesh
+from PySDM.physics import si
 
 
 def random_positions(grid: tuple, n_sd: int):
@@ -57,7 +58,7 @@ def test_strides(grid: tuple):
             Mesh(grid=(3, 4, 5, 6), size=(5, 5, 5, 5)),
         ),
     )
-    def test_cellural_attributes(mesh):
+    def test_cellular_attributes(mesh):
         # arrange
         n_sd = 666
         positions = random_positions(mesh.grid, n_sd)
@@ -87,3 +88,33 @@ def test_cellural_attributes(mesh):
                 0 <= min(position_in_cell[dim, :]) <= max(position_in_cell[dim, :]) < 1
             )
         assert position_in_cell.dtype == float
+
+    @staticmethod
+    @pytest.mark.parametrize(
+        "mesh, expected_dv",
+        (
+            (Mesh.mesh_0d(dv=44 * si.m**3), 44 * si.m**3),
+            (Mesh(grid=(100,), size=(1 * si.km,)), 10 * si.m * 1 * si.m * 1 * si.m),
+            (
+                Mesh(grid=(100, 200), size=(1 * si.km, 2 * si.km)),
+                10 * si.m * 10 * si.m * 1 * si.m,
+            ),
+            (
+                Mesh(grid=(10, 20, 30), size=(1 * si.km, 2 * si.km, 3 * si.km)),
+                (100 * si.m) ** 3,
+            ),
+        ),
+    )
+    def test_dv(mesh, expected_dv):
+        assert mesh.dv == expected_dv
+
+    @staticmethod
+    @pytest.mark.parametrize(
+        "mesh, expected_area",
+        (
+            (Mesh(grid=(100,), size=(44 * si.km,)), 1 * si.m**2),
+            (Mesh(grid=(100, 200), size=(44 * si.km, 666 * si.m)), 44 * si.km * si.m),
+        ),
+    )
+    def test_domain_bottom_surface_area(mesh, expected_area):
+        assert mesh.domain_bottom_surface_area == expected_area