refactor Paraview code & parcel-model animations (#1725)

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

Author: olastrz

.github/workflows/pdoc.yml

@@ -28,7 +28,7 @@ jobs:
           python-version: "3.12"
 
       - name: exclude pvpython scripts
-        run: rm examples/PySDM_examples/utils/pvanim.py
+        run: grep -rl '^#!/usr/bin/env pvpython' examples/PySDM_examples | xargs rm -f
 
       - run: pip install pdoc nbformat gitpython
       - run: pip install -e .

.github/workflows/tests.yml

@@ -300,8 +300,10 @@ jobs:
           rm /tmp/pytest/test_run_notebooks_*current
           ls /tmp/pytest/test_run_notebooks_*/
           ls /tmp/pytest/test_run_notebooks_*/output/
-          ffmpeg -i /tmp/pytest/test_run_notebooks_*/output/docs_intro_animation.ogv -loop 0 -vf "fps=10,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" $HOME/work/_temp/_github_home/figures/docs_intro_animation_${{ matrix.platform }}.gif
           mv /tmp/pytest/test_run_notebooks_*/output/last_animation_frame.pdf $HOME/work/_temp/_github_home/figures/last_animation_frame_${{ matrix.platform }}.pdf
+          for anim_name in docs_intro_animation parcel_animation; do
+              ffmpeg -i /tmp/pytest/test_run_notebooks_*/output/${anim_name}.ogv -loop 0 -vf "fps=10,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" $HOME/work/_temp/_github_home/figures/${anim_name}_${{ matrix.platform }}.gif
+          done
           ls $HOME/work/_temp/_github_home/figures
 
       - name: animation movie upload

PySDM/exporters/__init__.py

@@ -6,3 +6,4 @@
 from .netcdf_exporter_1d import NetCDFExporter_1d, readNetCDF_1d
 from .vtk_exporter import VTKExporter
 from .vtk_exporter_1d import VTKExporter_1d
+from .vtk_exporter_parcel import VTKExporterParcel

PySDM/exporters/vtk_exporter_parcel.py

@@ -0,0 +1,151 @@
+"""
+VTK Exporter for parcel PySDM simulations.
+
+This module defines `VTKExporterParcel`, a subclass of `PySDM.exporters.VTKExporter`,
+that writes simulation outputs to VTK format using `pyevtk`. It exports product
+profiles (e.g., relative humidity) as unstructured grids and particle attributes
+as point clouds, along with `.pvd` collection files for time-series visualization
+in ParaView.
+"""
+
+from pyevtk.hl import unstructuredGridToVTK, pointsToVTK
+from pyevtk.vtk import VtkHexahedron, VtkGroup
+import numpy as np
+
+from PySDM.exporters.vtk_exporter import VTKExporter
+
+
+# pylint: disable=too-many-instance-attributes
+class VTKExporterParcel(VTKExporter):
+    """
+    Custom VTK exporter for parcel PySDM, exporting products as grids
+    and attributes as point clouds for ParaView visualization.
+    """
+
+    def __init__(self, n_sd, output, mass_of_dry_air):
+        super().__init__()
+        self.output = output
+        self.coords = {
+            "x": np.random.random(n_sd),
+            "y": np.random.random(n_sd),
+            "z": np.random.random(n_sd),
+        }
+        self.half_diagonal = []
+        self.n_levels = len(self.output["products"]["z"])
+
+        volume = mass_of_dry_air / output["products"]["rhod"][0]
+        delta_z = output["products"]["z"][1] - output["products"]["z"][0]
+        for level in range(self.n_levels):
+            if level > 0:
+                prev_to_curr_density_ratio = (
+                    output["products"]["rhod"][level - 1]
+                    / output["products"]["rhod"][level]
+                )
+                volume *= prev_to_curr_density_ratio
+                delta_z = (
+                    output["products"]["z"][level] - output["products"]["z"][level - 1]
+                )
+            area = volume / delta_z
+            self.half_diagonal.append((2 * area) ** 0.5)
+
+    def write_pvd(self):
+        pvd_attributes = VtkGroup(self.attributes_file_path)
+        for key, value in self.exported_times["attributes"].items():
+            pvd_attributes.addFile(key + ".vtu", sim_time=value)
+        pvd_attributes.save()
+
+        pvd_products = VtkGroup(self.products_file_path)
+        for key, value in self.exported_times["products"].items():
+            pvd_products.addFile(key + ".vtu", sim_time=value)
+        pvd_products.save()
+
+    def export_products(
+        self, step, simulation
+    ):  # pylint: disable=arguments-differ, too-many-locals
+        path = self.products_file_path + "_num" + self.add_leading_zeros(step)
+        self.exported_times["products"][path] = step * simulation.particulator.dt
+
+        n_vertices = 4 * self.n_levels
+        x_vertices = np.zeros(n_vertices)
+        y_vertices = np.zeros(n_vertices)
+        z_vertices = np.zeros(n_vertices)
+        connectivity = [0, 1, 2, 3]
+        cell_type = np.zeros(self.n_levels - 1)
+        cell_type[:] = VtkHexahedron.tid
+
+        for level in range(self.n_levels):
+            half_diag = self.half_diagonal[level]
+            current_z = self.output["products"]["z"][level]
+            idx = level * 4
+            x_vertices[idx], y_vertices[idx], z_vertices[idx] = (
+                -half_diag,
+                -half_diag,
+                current_z,
+            )
+            idx += 1
+            x_vertices[idx], y_vertices[idx], z_vertices[idx] = (
+                -half_diag,
+                half_diag,
+                current_z,
+            )
+            idx += 1
+            x_vertices[idx], y_vertices[idx], z_vertices[idx] = (
+                half_diag,
+                half_diag,
+                current_z,
+            )
+            idx += 1
+            x_vertices[idx], y_vertices[idx], z_vertices[idx] = (
+                half_diag,
+                -half_diag,
+                current_z,
+            )
+            connectivity += [*range(4 * (level + 1), 4 * (level + 2))] * 2
+        connectivity = np.asarray(connectivity[:-4])
+        offset = np.asarray(range(8, 8 * self.n_levels, 8))
+
+        _ = {"test_pd": np.array([44] * n_vertices)}
+
+        relative_humidity = self.output["products"]["S_max_percent"]
+        cell_data = {
+            "RH": np.full(shape=(len(relative_humidity) - 1,), fill_value=np.nan)
+        }
+        cell_data["RH"][:step] = (
+            np.array(relative_humidity[:-1] + np.diff(relative_humidity) / 2)
+        )[:step]
+        unstructuredGridToVTK(
+            path,
+            x_vertices,
+            y_vertices,
+            z_vertices,
+            connectivity=connectivity,
+            offsets=offset,
+            cell_types=cell_type,
+            cellData=cell_data,
+        )
+
+    def export_attributes(self, step, simulation):  # pylint: disable=arguments-differ
+        path = self.attributes_file_path + "_num" + self.add_leading_zeros(step)
+        self.exported_times["attributes"][path] = step * simulation.particulator.dt
+
+        payload = {}
+
+        for key in self.output["attributes"].keys():
+            payload[key] = np.asarray(self.output["attributes"][key])[:, step].copy()
+        if step != 0:
+            delta_z = (
+                self.output["products"]["z"][step]
+                - self.output["products"]["z"][step - 1]
+            )
+            if step == 1:
+                self.coords["z"] *= delta_z
+            else:
+                self.coords["z"] += delta_z
+
+        pointsToVTK(
+            path,
+            2 * (self.coords["x"] - 0.5) * self.half_diagonal[step],
+            2 * (self.coords["y"] - 0.5) * self.half_diagonal[step],
+            self.coords["z"],
+            data=payload,
+        )

examples/PySDM_examples/Pyrcel/simulation.py

@@ -13,7 +13,13 @@
 
 class Simulation(BasicSimulation):
     def __init__(
-        self, settings, products=None, scipy_solver=False, rtol_thd=1e-10, rtol_x=1e-10
+        self,
+        settings,
+        products=None,
+        scipy_solver=False,
+        rtol_thd=1e-10,
+        rtol_x=1e-10,
+        mass_of_dry_air=44 * si.kg,
     ):
         n_sd = sum(settings.n_sd_per_mode)
         builder = Builder(
@@ -27,7 +33,7 @@ def __init__(
                 initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio,
                 T0=settings.initial_temperature,
                 w=settings.vertical_velocity,
-                mass_of_dry_air=44 * si.kg,
+                mass_of_dry_air=mass_of_dry_air,
             ),
         )
         builder.add_dynamic(AmbientThermodynamics())

examples/PySDM_examples/Strzabala_2025_BEng/__init__.py

@@ -0,0 +1,8 @@
+# pylint: disable=invalid-name
+"""
+ParaView visualisation of a parcel-model example based on the test case from
+[Pyrcel package docs](https://pyrcel.readthedocs.io/)
+
+paraview_parcel_model.ipynb:
+.. include:: ./paraview.ipynb.badges.md
+"""