VTK algorithms of the temporal filters, in situ and driving

changelog.d/vtk-algorithm-cache-observers.md

@@ -2,7 +2,10 @@
 
 - Python crashed in its garbage collector when it freed, in a reference cycle (a
   traceback's frames, say), a VTK dataset one of Tack's VTK algorithms had read: the
-  observer that drops the dataset's cache entry, referenced only by the dataset, was
-  cleared as garbage before the deletion called it (a VTK bug, in 9.7.1 and master).
-  The observers are kept reachable until they fire, and only note the deletion; the
-  cache drops the entry at its next use.
+  `DeleteEvent` observer that dropped the dataset's cache entry, referenced only by
+  the dataset, was cleared as garbage before the deletion called it (a VTK bug, in
+  9.7.1 and master). The cache entry is now an attribute of the dataset's Python
+  wrapper, which VTK keeps while the dataset lives and frees with it: no observer.
+  An entry is converted again when the dataset's modification time or time step
+  changes; after changing a dataset in place, call `Modified()`, as VTK's own
+  filters need.

changelog.d/vtk-algorithms-amr.md

@@ -0,0 +1,26 @@
+### Added
+
+- The VTK algorithms take an AMR hierarchy (`vtkOverlappingAMR`) whole, as one
+  `AMRDataSet`, where the Tack filter takes one: `Probe` (as its source),
+  `Statistics`, `Entropy`, `Integrate`, `Voxelize` and `GhostCells`, and the new
+  `AverageDown` (`td.average_down`, its input left as it was). A region that a
+  finer box covers, and a point that a sibling box holds, counts once.
+  `GhostCells` and `AverageDown` give a `vtkOverlappingAMR` of the input's boxes,
+  and the next AMR filter takes the hierarchy as it was made, ghost layers too.
+  The others run box by box, as Tack's filters do on AMR, each level's boxes as
+  pieces and the cells finer boxes cover left out: a field filter gives the
+  hierarchy with its fields, one that makes a dataset a
+  `vtkPartitionedDataSetCollection` with a partitioned dataset per level.
+  `GenerateGlobalIds`, `Streamlines` and `Advect` refuse a hierarchy. Over MPI
+  ranks each rank holds its own blocks, as VTK spreads a hierarchy.
+- `vtk_to_dataset(..., world=)`: the ranks a `vtkOverlappingAMR`'s blocks are
+  spread over, a `World` or an mpi4py communicator, by default VTK's (`world()`).
+
+### Fixed
+
+- A `vtkOverlappingAMR` given to a VTK algorithm gave an empty hierarchy, with no
+  error: VTK makes a hierarchy a partitioned collection, which they filtered as a
+  tree into an output of the hierarchy's type that could not hold it.
+- `vtk_to_dataset` of a `vtkOverlappingAMR` spread over MPI ranks marked each
+  rank's boxes as if its own were the whole hierarchy; it now converts in VTK's
+  world unless given one.

changelog.d/vtk-algorithms-temporal.md

@@ -0,0 +1,21 @@
+### Added
+
+- `RunningStatistics`, `TimeDerivative`, `Advection`, `Pathlines` and
+  `Streaklines`, VTK algorithms of Tack's temporal filters, in two modes. In
+  situ, each execution takes the step the input holds and gives the output at
+  it, stamped with its time; the same step again is not taken twice, and a step
+  back in time, a parameter changed or VTK's `NO_PRIOR_TEMPORAL_ACCESS` reset
+  starts over. Driving, the filter asks its input for the steps it needs, one
+  per pass and never one it has taken: the statistics, particles and paths over
+  every step (or those in `time_range`) in one output with no time (time in the
+  data: `step` and `time` on the particles, `time` on the paths' points), and
+  streaklines and derivatives at the time requested, going on from the steps
+  taken. A filter drives where its input has time steps and the pipeline can go
+  back in time; `mode=` chooses.
+
+### Fixed
+
+- `Advection` (and `Pathlines`, `Streaklines`) over pieces on MPI ranks waited
+  forever where a rank had no particles, released particles at other steps, or
+  saw its particles stop at another substep: the probes are collective, and a
+  rank skipped them. Every rank now takes part in each.

docs/data/amr-and-octrees.md

@@ -132,7 +132,8 @@ every rank, each as its `(lo, hi)` cell indices.
 back for images and rectilinear boxes. A hierarchy with curvilinear boxes
 goes to VTK as a `vtkPartitionedDataSetCollection` with one partitioned
 dataset per level. `td.amr_gaussian_pulse` makes VTK's Gaussian-pulse test
-hierarchy.
+hierarchy. In a VTK pipeline, Tack's AMR filters take a `vtkOverlappingAMR`
+whole (see [VTK](../visualization/vtk.md)).
 
 ## Tree-based AMR
 

docs/visualization/over-time.md

@@ -270,8 +270,63 @@ state, so every step must bring the same pieces, and the output is pieces
 too. Over MPI ranks each rank feeds its own pieces. The statistics and
 derivatives need nothing from other ranks. `Advection` over pieces probes
 each point on the rank whose cells hold it, so particles stay on the rank
-that released them. [Filters over pieces](over-pieces.md) explains pieces
-and ranks.
+that released them. Every rank takes part in each probe, so a rank may hold
+no seeds at all. [Filters over pieces](over-pieces.md) explains pieces and
+ranks.
+
+## In a VTK pipeline
+
+`tack.interop.vtk.algorithms` has each of these filters as a VTK algorithm:
+`RunningStatistics`, `TimeDerivative`, `Advection`, `Pathlines` and
+`Streaklines`. Each works in one of two modes.
+
+- **In situ**, the filter cannot drive the pipeline. Each execution takes the
+  step its input holds and gives the output at that step, stamped with its
+  time: the statistics or the paths up to that step, the particles at it. An
+  execution at the time of the step it took last gives the same output again
+  without taking the step twice. A step back in time, a changed parameter, or
+  VTK's reset signal starts the filter over.
+- **Driving**, the filter asks its input for the steps it needs, one per pass,
+  and never asks again for a step it has taken. `RunningStatistics`,
+  `Advection` and `Pathlines` take every step (or those in `time_range`) and
+  give one output after the last, with no time on it. The time is in the
+  data: the paths' `time` on their points, and for `Advection` the particles
+  at every step in one cloud, with `step` and `time` on the points.
+  `Streaklines` and `TimeDerivative` give their answer at the time requested
+  of them, stamped with it. They go on from the steps they have taken, and a
+  request for an earlier time starts them over.
+
+A filter drives where its input has time steps, unless VTK's
+`NO_PRIOR_TEMPORAL_ACCESS` key says the pipeline cannot go back in time, as
+an in situ adaptor sets it. `mode="in situ"` or `mode="driving"` chooses
+either. Driving to a step gives what the same steps give in situ.
+
+<!-- tack-check: over-time-vtk -->
+```python
+from vtkmodules.vtkFiltersGeneral import vtkTimeSourceExample
+
+import tack
+from tack.interop.vtk.algorithms import RunningStatistics
+
+tack.set_default_device("cpu")
+
+# Driving: the filter runs the source through its ten steps itself.
+source = vtkTimeSourceExample()
+stats = RunningStatistics(fields=["Point Value"])
+stats.SetInputConnection(source.GetOutputPort())
+stats.Update()
+over_all = stats.GetOutputDataObject(0)
+assert over_all.GetFieldData().GetArray("count").GetValue(0) == 10
+
+# In situ: each step is handed to the filter, which gives the statistics so far.
+live = RunningStatistics(fields=["Point Value"], mode="in situ")
+for step in range(3):
+    source.UpdateTimeStep(step / 9)
+    live.SetInputDataObject(0, source.GetOutputDataObject(0))
+    live.Update()
+so_far = live.GetOutputDataObject(0)
+assert so_far.GetFieldData().GetArray("count").GetValue(0) == 3
+```
 
 ## Details
 

docs/visualization/vtk.md

@@ -214,6 +214,9 @@ their cell data.
 | Pieces and ghosts | `GhostCells`, `GenerateGlobalIds`, `MergePieces`, `GhostPoints` | over every piece of a `vtkPartitionedDataSet` |
 | | `Split` | a `vtkPartitionedDataSet` |
 | | `RemoveGhosts` | the input's structure |
+| AMR hierarchies | `AverageDown` | `vtkOverlappingAMR` |
+| Over time | `RunningStatistics`, `TimeDerivative` | fields, over the steps |
+| | `Advection`, `Pathlines`, `Streaklines` (seeds on port 1, or `seeds=`) | `vtkUnstructuredGrid` |
 
 A `vtkPartitionedDataSet` input runs the filter on each piece and gives a
 `vtkPartitionedDataSet`. The algorithms whose Tack filter answers for the
@@ -301,6 +304,60 @@ centroid, each field's integral on it, and the measure as `Length`, `Area` or
 `Volume`; integrals are by quadrature through each cell's map, so curved and
 higher-order cells exactly.
 
+An AMR hierarchy, a `vtkOverlappingAMR`, is one Tack `AMRDataSet` (see
+[AMR and octrees](../data/amr-and-octrees.md)), each box marked with what a
+finer box covers. The algorithms whose Tack filter takes a hierarchy take it
+whole: `Probe` as its source (each point at its finest box, as
+`vtkAMRResampleFilter` samples), `Statistics`, `Entropy`, `Integrate`,
+`Voxelize`, `GhostCells` and `AverageDown`. A region that a finer box covers,
+and a point that a sibling box holds, counts once. `GhostCells` grows each box
+within its level and gives a `vtkOverlappingAMR` of the input's boxes whose
+blocks are the grown ones. `AverageDown` puts the finer levels' values of one
+field into the coarse cells and points they cover and leaves its input as it
+was.
+
+The other algorithms run box by box, as Tack's filters do on AMR, each level's
+boxes taken as the pieces of a partitioned dataset and the cells a finer box
+covers left out. One that computes fields (`Elevation`, `Gradients`,
+`Project`, ...) gives the hierarchy back with them. One that makes a dataset
+(`Contour`, `Clip`, `Slice`, `Threshold`, ...) gives a
+`vtkPartitionedDataSetCollection` with a partitioned dataset per level, named
+`level0`, `level1`, and so on. Contours and slices of different levels do not
+join where the levels meet ([Adaptive data](adaptive-data.md) has the limits),
+and `Gradients` at a box's edge needs `GhostCells` first, as for any pieces.
+`GenerateGlobalIds`, `Streamlines` and `Advect` refuse a hierarchy, since box by
+box their answers would mislead: ids repeated from level to level, paths
+ending where a finer box begins. Over MPI ranks, each rank holds its own blocks, the others'
+empty, as VTK spreads a hierarchy, and the ranks convert it together
+(`vtk_to_dataset(amr, world=...)` does the same).
+
+<!-- tack-check: vtk-pipeline-amr -->
+```python
+from vtkmodules.vtkFiltersAMR import vtkAMRGaussianPulseSource
+
+import tack
+from tack.interop.vtk.algorithms import AverageDown, GhostCells, Statistics, Threshold
+
+tack.set_default_device("cpu")
+
+pulse = vtkAMRGaussianPulseSource(dimension=3)()       # a box, two finer ones in it
+once = Statistics(fields=["Gaussian-Pulse"])(pulse)    # covered cells counted once
+
+grown = GhostCells(layers=1)(pulse)                    # blocks grown by a layer
+assert grown.GetDataSetAsCartesianGrid(1, 0).number_of_cells > 4 * 6 * 10
+again = Statistics(fields=["Gaussian-Pulse"])(grown)   # ghosts count nowhere
+assert again.GetColumnByName("n").GetValue(0) == once.GetColumnByName("n").GetValue(0)
+
+averaged = AverageDown(field="Gaussian-Pulse")(pulse)  # fine values in the coarse box
+
+kept = Threshold(field=("Gaussian-Pulse", "cells"), lower=0.0, upper=1.0)(pulse)
+assert kept.GetNumberOfPartitionedDataSets() == 2       # a level each, box by box
+```
+
+The temporal filters run in situ, one step per execution, or drive the
+pipeline through the input's time steps themselves. See
+[Steps over time](over-time.md#in-a-vtk-pipeline).
+
 ## Write your own VTK algorithm
 
 Subclass `TackAlgorithm` and give it `compute(data)`, a Tack dataset from a

packages/tack-vis/src/tack/dataset/temporal.py

@@ -372,12 +372,15 @@ def _frame_of(data):
 
 def _velocities(data, field, x):
     """The velocity at world points ``x`` (``(n, 3)`` f64) and whether a cell holds
-    each: ``probe``, over pieces routed to the rank whose cells hold them."""
+    each: ``probe``, over pieces routed to the rank whose cells hold them. Over
+    pieces the probe is collective, so a rank with no points probes too."""
     from tack.dataset.algorithms.probe import probe
 
-    if not len(x):
+    if not len(x) and not isinstance(data, PartitionedDataSet):
         return np.zeros((0, 3)), np.zeros(0, bool)
     out = probe(data, x, fields=[field])
+    if not len(x):
+        return np.zeros((0, 3)), np.zeros(0, bool)
     values = np.asarray(out.fields[field].values.numpy(), np.float64).reshape(len(x), -1)
     return values, out.fields["valid"].values.numpy() == 1
 
@@ -471,7 +474,7 @@ def _velocity_scale(self, data):
 
     def _release(self, data, t1):
         """Particles pending with a release time at or before ``t1``."""
-        keep = []
+        keep, due = [], []
         for points, time, sources in self._pending:
             time = (self._last if self._last is not None else t1) if time is None else time
             if time > t1:
@@ -480,15 +483,30 @@ def _release(self, data, t1):
             if self._last is None and time < t1:
                 raise ValueError(f"particles released at {time}, before the first step's "
                                  f"time {t1}")
+            due.append((points, time, sources))
+        # Released outside any cell: SUCCESS cleared, as advect's seeds. Those
+        # released before t1 are probed on the previous step, the rest on this one,
+        # each group at once: over pieces every rank probes both, its points or none,
+        # since the probe is collective (a rank with fewer releases waited forever).
+        earlier = [self._previous is not None and time < t1 for _, time, _ in due]
+        inside = [None] * len(due)
+        for where, which in ((self._previous, True), (data, False)):
+            group = [i for i, e in enumerate(earlier) if e == which]
+            if where is None or not (group or isinstance(data, PartitionedDataSet)):
+                continue
+            points = np.concatenate([due[i][0] for i in group]) if group else \
+                np.zeros((0, 3))
+            _, found = _velocities(where, self.field, points)
+            for i, part in zip(group, np.split(found, np.cumsum(
+                    [len(due[i][0]) for i in group])[:-1]) if group else []):
+                inside[i] = part
+        for (points, time, sources), found in zip(due, inside):
             n = len(points)
             ids = self._released + np.arange(n, dtype=np.int64)
             if isinstance(data, PartitionedDataSet):
                 ids += np.int64(data.world.rank) << np.int64(40)
             self._released += n
-            # Released outside any cell: SUCCESS cleared, as advect's seeds.
-            where = self._previous if self._previous is not None and time < t1 else data
-            _, inside = _velocities(where, self.field, points)
-            status = np.where(inside, 1, 0).astype(np.int32)
+            status = np.where(found, 1, 0).astype(np.int32)
             for name, values in (("x", points), ("time", np.full(n, time)),
                                  ("release", np.full(n, time)), ("seed", ids),
                                  ("source", sources), ("status", status)):
@@ -518,10 +536,14 @@ def velocity(x, t):
             return (1.0 - a) * v0 + a * v1, in0 & in1
 
         limit = self.step_size if self.step_size is not None else span
+        # Over pieces every rank takes each substep, moving particles or not, until
+        # none moves on any: the probes are collective.
+        world = now.world if isinstance(now, PartitionedDataSet) else None
         while True:
             moving = np.flatnonzero(((self.status & 1) != 0) & ((self.status & 2) == 0)
                                     & (self.time < t1 - 1e-12 * max(abs(t1), 1.0)))
-            if not len(moving):
+            if not (len(moving) if world is None
+                    else any(world._allgather(bool(len(moving))))):
                 break
             x, t = self.x[moving], self.time[moving]
             h = np.minimum(limit, t1 - t)

packages/tack-vis/src/tack/interop/vtk/__init__.py

@@ -986,7 +986,7 @@ def __repr__(self):
 
 
 def vtk_to_dataset(grid, dtype=_DEFAULT_FLOAT, polyhedral=None, orient=False, id_dtype=None,
-                   origin=None, copy=True, merge_points=None):
+                   origin=None, copy=True, merge_points=None, world=None):
     """Copy a ``vtkUnstructuredGrid``, ``vtkPolyData``, ``vtkStructuredGrid``,
     ``vtkRectilinearGrid`` or ``vtkImageData`` into a ``tack.dataset.DataSet``: an image's points as
     ``UniformCoordinates`` (with its direction matrix), a rectilinear grid's as a
@@ -1056,7 +1056,10 @@ def vtk_to_dataset(grid, dtype=_DEFAULT_FLOAT, polyhedral=None, orient=False, id
     rectilinear grid) converted as above, placed in its level by its
     ``vtkAMRBox``, the levels' refinement ratios kept; the refined cells are
     marked afresh (VTK's ``vtkGhostType`` bits come in and are recomputed).
-    Blocks with ghost layers are refused.
+    Blocks with ghost layers are refused. Over MPI ranks each holds its own
+    blocks (the others' empty, as VTK spreads a hierarchy) and every rank
+    converts together: ``world`` is the ranks, a ``tack.communication.World`` or
+    an mpi4py communicator, by default VTK's (``world()``).
 
     A ``vtkCellGrid`` (VTK's discontinuous-Galerkin cells, read from ``.dg``
     files by ``vtkCellGridReader``) becomes a ``DataSet`` of its cells on its
@@ -1114,7 +1117,7 @@ def vtk_to_dataset(grid, dtype=_DEFAULT_FLOAT, polyhedral=None, orient=False, id
         origin = _far_origin(grid)
     out = _vtk_to_dataset(grid, dtype=dtype, polyhedral=polyhedral, orient=orient,
                           id_dtype=id_dtype, origin=origin, copy=copy,
-                          merge_points=merge_points)
+                          merge_points=merge_points, world=world)
     if grid.GetInformation().Has(_time_key()):
         _stamp(out, grid.GetInformation().Get(_time_key()))
     return out
@@ -1407,7 +1410,8 @@ def _take_cell_data(source, target, rows):
 
 
 def _vtk_to_dataset(grid, dtype=tack.f32, polyhedral=None, orient=False, id_dtype=None,
-                   origin=None, copy=True, merge_points=None, placement=None, _device=None):
+                   origin=None, copy=True, merge_points=None, placement=None, _device=None,
+                   world=None):
     """``vtk_to_dataset`` before the time is stamped. ``placement`` is a structured
     grid's offset in a larger grid (three values), in place of its extent's start:
     an AMR box's place in its level. ``_device`` is the device a ``copy=False``
@@ -1444,7 +1448,7 @@ def _vtk_to_dataset(grid, dtype=tack.f32, polyhedral=None, orient=False, id_dtyp
 
     shift = np.asarray(as_origin(origin))
     if isinstance(grid, vtkOverlappingAMR):
-        return _amr_to_dataset(grid, options)
+        return _amr_to_dataset(grid, options, world)
     from vtkmodules.vtkCommonDataModel import vtkPolyData
 
     if isinstance(grid, vtkPolyData):
@@ -1732,10 +1736,18 @@ def _to_vtk_ids(values):
     return np.ascontiguousarray(np.asarray(values), np.int64)
 
 
-def _amr_to_dataset(amr, options):
+def _amr_to_dataset(amr, options, world=None):
     """A ``vtkOverlappingAMR`` as an ``AMRDataSet``: each level's blocks, each placed
-    in its level by its ``vtkAMRBox``."""
+    in its level by its ``vtkAMRBox``; this rank's, in ``world`` (``vtk_to_dataset``'s,
+    VTK's by default)."""
+    from tack.communication import World
     from tack.dataset import AMRDataSet, StructuredTopology
+    from tack.interop.vtk import world as vtk_world
+
+    if world is None:
+        world = vtk_world()
+    elif not isinstance(world, World):
+        world = World(world)                # an mpi4py communicator
 
     levels, ratios = [], []
     for level in range(amr.GetNumberOfLevels()):
@@ -1764,7 +1776,7 @@ def _amr_to_dataset(amr, options):
         levels.append(boxes)
         if level + 1 < amr.GetNumberOfLevels():
             ratios.append(int(amr.GetRefinementRatio(level)))
-    return AMRDataSet(levels, ratios)
+    return AMRDataSet(levels, ratios, world=world)
 
 
 def _bits(array):