diff --git a/changelog.d/vtk-algorithm-cache-observers.md b/changelog.d/vtk-algorithm-cache-observers.md index 9edde338a..37e86f3ab 100644 --- a/changelog.d/vtk-algorithm-cache-observers.md +++ b/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. diff --git a/changelog.d/vtk-algorithms-amr.md b/changelog.d/vtk-algorithms-amr.md new file mode 100644 index 000000000..8501d507a --- /dev/null +++ b/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. diff --git a/changelog.d/vtk-algorithms-temporal.md b/changelog.d/vtk-algorithms-temporal.md new file mode 100644 index 000000000..31ec823e4 --- /dev/null +++ b/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. diff --git a/docs/data/amr-and-octrees.md b/docs/data/amr-and-octrees.md index b34411409..ab7f05eea 100644 --- a/docs/data/amr-and-octrees.md +++ b/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 diff --git a/docs/visualization/over-time.md b/docs/visualization/over-time.md index b263c05dd..f9176e6d6 100644 --- a/docs/visualization/over-time.md +++ b/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. + + +```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 diff --git a/docs/visualization/vtk.md b/docs/visualization/vtk.md index 9086b97a6..78cd504e5 100644 --- a/docs/visualization/vtk.md +++ b/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). + + +```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 diff --git a/packages/tack-vis/src/tack/dataset/temporal.py b/packages/tack-vis/src/tack/dataset/temporal.py index 0aa98d2a3..bf5dee2e5 100644 --- a/packages/tack-vis/src/tack/dataset/temporal.py +++ b/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) diff --git a/packages/tack-vis/src/tack/interop/vtk/__init__.py b/packages/tack-vis/src/tack/interop/vtk/__init__.py index 0f3e0fa54..9d26c82dc 100644 --- a/packages/tack-vis/src/tack/interop/vtk/__init__.py +++ b/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): diff --git a/packages/tack-vis/src/tack/interop/vtk/algorithms.py b/packages/tack-vis/src/tack/interop/vtk/algorithms.py index c5c05e845..ce803580f 100644 --- a/packages/tack-vis/src/tack/interop/vtk/algorithms.py +++ b/packages/tack-vis/src/tack/interop/vtk/algorithms.py @@ -28,6 +28,19 @@ block of its own, and the output is a tree of the same shape, with the blocks' names and a collection's data assembly. +An AMR hierarchy (``vtkOverlappingAMR``) is one Tack ``AMRDataSet``, its boxes marked +with what a finer box covers (``REFINEDCELL``) and what a sibling holds. The filters +whose Tack algorithm takes one take it whole -- ``Probe`` (as its source), +``Statistics``, ``Entropy``, ``Integrate``, ``Voxelize``, ``GhostCells`` and +``AverageDown`` -- so each region counts once; a hierarchy out (``GhostCells``, +``AverageDown``) has the input's boxes. The others run box by box, each level's +boxes as the pieces of a partitioned dataset, leaving out the cells finer boxes +cover, as Tack's filters do on AMR: a filter that computes fields gives the +hierarchy with them, one that makes a dataset a ``vtkPartitionedDataSetCollection`` +with a partitioned dataset per level (``level0``, ...), whose contours and slices do +not join where the levels meet. ``GenerateGlobalIds``, ``Streamlines`` and +``Advect``, whose answer box by box would mislead, refuse one. + Filters that make a dataset give a ``vtkUnstructuredGrid`` (or keep a structured grid's type where the Tack filter does); filters that compute a field (``Gradients``, ``Elevation``, ``CellQuality``) give their input's type, shallow-copied, with the @@ -46,6 +59,11 @@ function, so the same zero set and values; a ``vtkPlane`` is normalized, its values those of the unit normal. +The temporal filters (``RunningStatistics``, ``TimeDerivative``, ``Advection``, +``Pathlines``, ``Streaklines``) run in situ, taking the step their input holds at +each execution and giving the output at it, or drive the pipeline through their +input's time steps, one per pass and none twice: ``_Temporal``. + An exception in a filter is VTK's error: printed, the output left empty and ``Update()`` false. It is also kept as the algorithm's ``error`` (``None`` after a run that succeeded), for Python callers. @@ -55,8 +73,10 @@ __all__ = [ "Advect", + "Advection", "AsPolygons", "AsPolyhedra", + "AverageDown", "BoundaryFaces", "CellCenters", "CellGeometry", @@ -102,11 +122,13 @@ "MergePieces", "MergePoints", "Outliers", + "Pathlines", "PointIds", "PointNormals", "Probe", "Project", "RemoveGhosts", + "RunningStatistics", "Shrink", "Slice", "SlicePlane", @@ -115,6 +137,7 @@ "Spherical", "Split", "Statistics", + "Streaklines", "Streamlines", "Subdivide", "SurfaceNormals", @@ -122,6 +145,7 @@ "Tetrahedralize", "Threshold", "ThresholdPoints", + "TimeDerivative", "Transform", "Triangulate", "Tubes", @@ -139,8 +163,10 @@ vtkDataObject, vtkMultiBlockDataSet, vtkMultiPieceDataSet, + vtkOverlappingAMR, vtkPartitionedDataSet, vtkPartitionedDataSetCollection, + vtkUniformGridAMR, ) from vtkmodules.vtkCommonExecutionModel import vtkAlgorithm @@ -154,72 +180,120 @@ # A Tack dataset keeps what it derives from its topology -- locators, faces, point # links, groups -- on the topology, so filters that read one VTK dataset (several # on one input, or one run again after a parameter changed) build that once when -# they share one Tack dataset. Kept by the VTK object's address, modification time -# and time step: the modification time is a global counter, so the same object -# modified, or a later one at the same address, does not match. An entry goes when -# its VTK object is deleted, and holds nothing that keeps that object alive. Beside -# the dataset sharing VTK's memory, an entry keeps its copy on each device a filter -# asked for (``device=``), so the data moves there once. +# they share one Tack dataset. The entry is an attribute of the VTK dataset's +# Python wrapper (``_CACHE``): VTK keeps a wrapper's attributes while its C++ object +# lives, on every later wrapper of it (the one GetInputData gives in the next +# filter), and frees them with the object, so an entry lives exactly as long as its +# dataset, with no observer to call (a DeleteEvent observer crashed Python's +# collector, VTK bugs 19 in the vault). On VTK master a wrapper costs the same however +# many objects keep attributes (0.6 us); VTK 9.7.1 sweeps them at every wrapper (50 us +# with 10,000), as it does the attributes the pythonic API gives every wrapper. +# +# An entry is stale when the dataset's modification time or time step differs from +# what it was cached at (``_key``). The modification time follows VTK's contract: it +# sees a change made with ``Modified()`` (an array's or the points'), arrays added or +# removed, not values written in place without it, nor ``InsertNextCell``. Tack shares +# VTK's memory, so such values reach its fields, but what it derived from the old +# state (a locator, the faces) is kept: call ``Modified()`` after changing a dataset +# in place, as VTK's own filters need. Beside the dataset sharing VTK's memory, an +# entry keeps its copy on each device a filter asked for (``device=``), so the data +# moves there once. +# +# An AMR hierarchy's entry is kept on its metadata (``_holder``), which VTK shares +# between a hierarchy and its shallow copies -- what calling a filter gives is one, +# made of new blocks over the same arrays -- and is valid for its blocks' arrays, +# each by address and modification time (``_signature``): the hierarchy's own +# modification time does not follow its blocks'. A hierarchy ``GhostCells`` grew +# reaches the next filter so, which could not convert it again. -_datasets = {} # address -> ((modification time, time step), {device: DataSet}) +#: The attribute of a VTK dataset's wrapper that holds its entry: ``((modification +#: time, time step), {device: DataSet})``. +_CACHE = "_tack_cache" -#: Addresses of VTK datasets deleted since the cache was last used. A DeleteEvent only -#: notes one, and its entry goes at the next use (``_drop_deleted``), outside whatever -#: freed the dataset. -_deleted = [] -#: Each DeleteEvent observer, by its dataset's address, kept here until it has fired. -#: Referenced only by its VTK object, an observer is garbage when Python's collector -#: frees that object in a cycle (a traceback's frame holding it): the collector clears -#: the function's globals and closure, then the deletion calls it, and Python crashed. -_observers = {} +def _key(vtk_data): + """What an entry is valid for: the dataset's modification time and time step; an + AMR hierarchy's blocks' arrays (``_signature``) and time step.""" + from tack.interop.vtk import _time_key + info = vtk_data.GetInformation() + time = info.Get(_time_key()) if info.Has(_time_key()) else None + if isinstance(vtk_data, vtkUniformGridAMR): + return _signature(vtk_data), time + return vtk_data.GetMTime(), time + + +def _holder(vtk_data): + """What keeps ``vtk_data``'s entry: its wrapper; an AMR hierarchy's metadata, which + its shallow copies share.""" + if isinstance(vtk_data, vtkOverlappingAMR): + meta = vtk_data.GetOverlappingAMRMetaData() + if meta is not None: + return meta + return vtk_data + + +def _signature(amr): + """An AMR hierarchy as its entry knows it: each block here, its place and lattice + (an image's origin and spacing) and its arrays (a rectilinear grid's coordinates + too) by address and modification time. A shallow copy has the same; a block + replaced, an array changed with ``Modified()``, added or removed, does not. The + modification time is VTK's global counter, so a new array at a freed one's + address differs.""" + from vtkmodules.vtkCommonDataModel import vtkRectilinearGrid + + def stamp(array): + return array.GetName(), array.GetAddressAsString("vtkObject"), array.GetMTime() -def _note_deleted(address): - """A DeleteEvent observer noting ``address`` deleted.""" - def note(caller, event): - _deleted.append(address) - return note + out = [] + for level, index, block in _blocks(amr): + if isinstance(block, vtkRectilinearGrid): + lattice = tuple(stamp(a) for a in (block.GetXCoordinates(), + block.GetYCoordinates(), + block.GetZCoordinates())) + else: + lattice = block.GetOrigin(), block.GetSpacing() + arrays = tuple(stamp(data.GetAbstractArray(i)) + for data in (block.GetPointData(), block.GetCellData()) + for i in range(data.GetNumberOfArrays())) + out.append((level, index, block.GetDimensions(), lattice, arrays)) + return tuple(out) -def _drop_deleted(): - """Drop the entries, and observers, of the VTK datasets deleted since the last use.""" - while _deleted: - address = _deleted.pop() - _datasets.pop(address, None) - _observers.pop(address, None) +def _blocks(amr): + """``(level, index, block)`` for each block an AMR hierarchy holds here (another + rank's are empty).""" + for level in range(amr.GetNumberOfLevels()): + for index in range(amr.GetNumberOfBlocks(level)): + block = amr.GetDataSetAsCartesianGrid(level, index) + if block is not None: + yield level, index, block -def _key(vtk_data): - from tack.interop.vtk import _time_key +def _remember(vtk_data, data): + """Keep ``data`` as the Tack dataset of ``vtk_data`` as it is now.""" + setattr(_holder(vtk_data), _CACHE, (_key(vtk_data), {None: data})) - info = vtk_data.GetInformation() - return vtk_data.GetMTime(), (info.Get(_time_key()) if info.Has(_time_key()) else None) +def _entry(vtk_data): + """``vtk_data``'s entry, or ``None`` (none, or a stale one).""" + entry = getattr(_holder(vtk_data), _CACHE, None) + return entry if entry is not None and entry[0] == _key(vtk_data) else None -def _remember(vtk_data, data): - """Keep ``data`` as the Tack dataset of ``vtk_data`` as it is now.""" - _drop_deleted() - address = vtk_data.GetAddressAsString("vtkObject") - if address not in _datasets: - # The address, not the caller: a DeleteEvent's caller is None. Only noted - # (see ``_deleted``); a later object at this address gets its own observer, - # since the noted address is dropped before it is remembered. - _observers[address] = observer = _note_deleted(address) - vtk_data.AddObserver("DeleteEvent", observer) - _datasets[address] = (_key(vtk_data), {None: data}) - - -def _as_dataset(vtk_data, device=None): + +def _as_dataset(vtk_data, device=None, world=None): """``vtk_to_dataset(vtk_data, dtype=None, copy=False)``, the same Tack dataset for the same VTK dataset unchanged (each piece of a ``vtkPartitionedDataSet``); - with ``device``, that dataset moved there (the same moved one each time).""" + with ``device``, that dataset moved there (the same moved one each time). An + AMR hierarchy is converted in ``world`` (``_as_hierarchy``).""" import tack from tack.dataset import PartitionedDataSet from tack.interop.vtk import _time_key, vtk_to_dataset if isinstance(device, str): device = tack.device(device) + if isinstance(vtk_data, vtkOverlappingAMR): + return _as_hierarchy(vtk_data, device, world) if isinstance(vtk_data, vtkPartitionedDataSet): pieces = [vtk_data.GetPartition(i) for i in range(vtk_data.number_of_partitions)] data = PartitionedDataSet(_as_dataset(piece, device) for piece in pieces @@ -230,11 +304,10 @@ def _as_dataset(vtk_data, device=None): data = vtk_to_dataset(vtk_data, dtype=None, copy=False) return data if device is None else data.to_device(device) return data - _drop_deleted() - entry = _datasets.get(vtk_data.GetAddressAsString("vtkObject")) - if entry is None or entry[0] != _key(vtk_data): + entry = _entry(vtk_data) + if entry is None: _remember(vtk_data, vtk_to_dataset(vtk_data, dtype=None, copy=False)) - entry = _datasets[vtk_data.GetAddressAsString("vtkObject")] + entry = getattr(vtk_data, _CACHE) copies = entry[1] data = copies[None] if device is None or data.device == device: @@ -244,6 +317,93 @@ def _as_dataset(vtk_data, device=None): return copies[device] +def _as_hierarchy(vtk_amr, device, world): + """A ``vtkOverlappingAMR`` as its ``AMRDataSet``, this rank's boxes sharing VTK's + memory, the hierarchy in ``world`` (by default VTK's): the same one while the + hierarchy and its blocks are unchanged, in the same world; with ``device``, its + boxes moved there.""" + from tack.interop.vtk import vtk_to_dataset + + entry = _entry(vtk_amr) + kept = entry is not None and (world is None or _same_world(entry[1][None].world, world)) + if world is not None and world.size > 1: + # Converting is collective: every rank converts again if one does. + kept = all(world._allgather(kept)) + if not kept: + _remember(vtk_amr, vtk_to_dataset(vtk_amr, dtype=None, copy=False, world=world)) + entry = _entry(vtk_amr) + copies = entry[1] + if device is None: + return copies[None] + if device not in copies: + copies[device] = _hierarchy_on(copies[None], device) + return copies[device] + + +def _same_world(a, b): + """Whether two Worlds are one: the same, or each this process alone (``world()`` + makes one each call without an MPI controller).""" + return a is b or (a.size == 1 and b.size == 1) + + +def _hierarchy_on(amr, device): + """``amr`` with its boxes on ``device``: itself where they are.""" + from tack.dataset import AMRDataSet + + levels = [[box.to_device(device) for box in pieces] for pieces in amr.levels] + if all(moved is box for moved, (_, box) in zip((b for lv in levels for b in lv), amr)): + return amr # memory the device addresses already + if amr.layers: + raise ValueError("an AMR hierarchy with ghost layers stays on its device: give " + "GhostCells the device instead") + return AMRDataSet._grown(amr, levels, amr._extents, amr.layers, amr._fills) + + +def _device_of(data): + """The device a Tack dataset, or an AMR hierarchy's first box, is on.""" + from tack.dataset import AMRDataSet + + if isinstance(data, AMRDataSet): + return next((box.device for _, box in data), None) + return data.device + + +def _amr_like(source): + """An empty ``vtkOverlappingAMR`` of ``source``'s hierarchy, with metadata of its own: + ``CopyStructure`` shares the input's, which holds the input's entry.""" + from vtkmodules.vtkCommonDataModel import vtkOverlappingAMRMetaData + + out = vtkOverlappingAMR() + out.CopyStructure(source) + meta = vtkOverlappingAMRMetaData() + meta.DeepCopy(source.GetOverlappingAMRMetaData()) + out.SetAMRMetaData(meta) + return out + + +def _amr_output(source, amr): + """A ``vtkOverlappingAMR`` of ``source``'s hierarchy (its boxes, origin, spacings and + ratios, on every rank) holding ``amr``'s boxes, sharing Tack's memory: this + rank's blocks of ``source`` in order, as they were converted. ``amr`` is + remembered as its Tack hierarchy: one ``GhostCells`` grew could not be converted + again, VTK's blocks with ghost layers being refused.""" + out = _amr_like(source) + held = {} + for level, index, _ in _blocks(source): + held.setdefault(level, []).append(index) + from tack.interop.vtk import dataset_to_vtk + + for level, pieces in enumerate(amr.levels): + indices = held.get(level, []) + if len(indices) != len(pieces): + raise ValueError(f"level {level} has {len(pieces)} boxes here, the input " + f"{len(indices)}") + for index, box in zip(indices, pieces): + out.SetDataSet(level, index, dataset_to_vtk(box, copy=False)) + _remember(out, amr) + return out + + class Parameter: """A property of an algorithm: read as the value set (or the default), set with ``Modified()`` so the pipeline runs again.""" @@ -263,6 +423,11 @@ def __get__(self, algorithm, owner=None): def __set__(self, algorithm, value): algorithm.__dict__[self.slot] = value + # Counted apart from the modification time, which the pipeline also moves (a + # new input given with SetInputDataObject): a temporal filter starts over + # when a parameter changes. + algorithm.__dict__["_parameter_changes"] = algorithm.__dict__.get( + "_parameter_changes", 0) + 1 algorithm.Modified() @@ -410,11 +575,16 @@ class TackAlgorithm(VTKPythonAlgorithmBase): filter that keeps the cells' structure); ``"points"``, a structured grid for any of those three and otherwise an unstructured grid (a filter that moves the points). ``INPUTS`` counts the input ports, the first ``REQUIRED_INPUTS`` of - them required.""" + them required. ``AMR_PORTS`` are the ports that take an AMR hierarchy + (``vtkOverlappingAMR``), the Tack function given the ``AMRDataSet`` whole; on + port 0 of any other filter it runs box by box (``_by_box``) unless ``BY_BOX`` + is false, and anywhere else one is refused (``_check_amr``).""" OUTPUT = "vtkUnstructuredGrid" INPUTS = 1 REQUIRED_INPUTS = None + AMR_PORTS = () + BY_BOX = True FUNCTION = None FIELDS = () IMPLICIT = () @@ -537,8 +707,8 @@ def _world(self): return kept[1] def _dataset(self, vtk_data, device): - """``_as_dataset``, a partitioned one in this run's World.""" - return _in_world(_as_dataset(vtk_data, device), self._run_world) + """``_as_dataset``, a partitioned one or a hierarchy in this run's World.""" + return _in_world(_as_dataset(vtk_data, device, self._run_world), self._run_world) def _supplied(self): """Arguments this run supplies where their parameters are unset: seeds read from @@ -554,16 +724,45 @@ def FillInputPortInformation(self, port, info): info.Set(vtkAlgorithm.INPUT_REQUIRED_DATA_TYPE(), _INPUT_TYPES[0]) for kind in _INPUT_TYPES[1:]: info.Append(vtkAlgorithm.INPUT_REQUIRED_DATA_TYPE(), kind) + if port in self.AMR_PORTS: + # A partitioned collection already, in VTK's model: named for what it is. + info.Append(vtkAlgorithm.INPUT_REQUIRED_DATA_TYPE(), "vtkOverlappingAMR") required = self.INPUTS if self.REQUIRED_INPUTS is None else self.REQUIRED_INPUTS if port >= required: info.Set(vtkAlgorithm.INPUT_IS_OPTIONAL(), 1) return 1 + def _check_amr(self, data, port): + """Refuse an AMR hierarchy where the filter takes none, whole or box by box. + VTK makes one a partitioned collection, which went through ``_tree`` level + by level into an output of the hierarchy's type that could not hold it: an + empty hierarchy, and no error.""" + if not isinstance(data, vtkUniformGridAMR): + return + if not isinstance(data, vtkOverlappingAMR): + raise TypeError(f"{type(self).__name__}: Tack's AMR is overlapping " + f"(vtkOverlappingAMR), not a {data.GetClassName()}") + if port not in self.AMR_PORTS and not (port == 0 and self.BY_BOX): + raise TypeError(f"{type(self).__name__} takes no AMR hierarchy on port {port}" + + ("" if port else ": box by box, its answer would mislead")) + + def _amr_output_type(self): + """The output's VTK class for a hierarchy filtered box by box.""" + return vtkPartitionedDataSetCollection + + def _check_inputs(self, inInfo): + for port in range(self.GetNumberOfInputPorts()): + for i in range(inInfo[port].GetNumberOfInformationObjects()): + self._check_amr(self.GetInputData(inInfo, port, i), port) + def _output_type(self, input): """The VTK class of the output for ``input``, a dataset or a partitioned dataset (trees keep their type).""" import vtkmodules.vtkCommonDataModel as dm + if isinstance(input, vtkOverlappingAMR): + # A hierarchy, unless the filter makes a type of its own. + return vtkOverlappingAMR if self.OUTPUT is None else getattr(dm, self.OUTPUT) if isinstance(input, vtkPartitionedDataSet): return vtkPartitionedDataSet if self.OUTPUT is None: @@ -583,7 +782,13 @@ def RequestDataObject(self, request, inInfo, outInfo): input = self.GetInputData(inInfo, 0, 0) if input is None: return 0 - kind = type(input) if isinstance(input, _TREES) else self._output_type(input) + if isinstance(input, vtkOverlappingAMR) and 0 in self.AMR_PORTS: + kind = self._output_type(input) + elif isinstance(input, vtkOverlappingAMR) and self.BY_BOX: + kind = self._amr_output_type() + else: + # A tree's type; a hierarchy this filter refuses (RequestData) its own. + kind = type(input) if isinstance(input, _TREES) else self._output_type(input) output = outInfo.GetInformationObject(0) current = output.Get(vtkDataObject.DATA_OBJECT()) if current is None or type(current) is not kind: @@ -611,20 +816,48 @@ def ProcessRequest(self, request, inInfo, outInfo): def RequestData(self, request, inInfo, outInfo): self._run_world = self._world() + self._check_inputs(inInfo) self._prepare(inInfo) input = self.GetInputData(inInfo, 0, 0) output = self.GetOutputData(outInfo, 0) - if isinstance(input, _TREES): + if isinstance(input, vtkOverlappingAMR) and 0 not in self.AMR_PORTS: + self._set_unit_output(output, self._by_box(input)) + elif isinstance(input, _TREES) and not isinstance(input, vtkOverlappingAMR): # Shares the blocks: what was made, and remembered, is the output's. output.CompositeShallowCopy(self._tree(input)) else: self._set_unit_output(output, self._unit(input)) return 1 + def _by_box(self, vtk_amr): + """The output for an AMR hierarchy this filter takes box by box: each level's + boxes, this rank's, as the pieces of a partitioned dataset (``_level``), the + cells finer boxes cover left out by the filter as their marks say; a + ``vtkPartitionedDataSetCollection`` of one partitioned dataset per level, + named ``level``.""" + from vtkmodules.vtkCommonDataModel import vtkCompositeDataSet + + from tack.interop.vtk import dataset_to_vtk + + amr = self._dataset(vtk_amr, self.device) + out = vtkPartitionedDataSetCollection() + out.SetNumberOfPartitionedDataSets(amr.num_levels) + for level, pieces in enumerate(amr.levels): + result = self._level(pieces) + out.SetPartitionedDataSet(level, vtkPartitionedDataSet() if result is None else + _as_partitioned(dataset_to_vtk(result, copy=False))) + out.GetMetaData(level).Set(vtkCompositeDataSet.NAME(), f"level{level}") + return out + + def _level(self, pieces): + """The result for one level's boxes, a ``PartitionedDataSet``.""" + return self.execute(pieces) + def _prepare(self, inInfo): """Read what every unit of a run needs from the other ports (a source).""" def _set_unit_output(self, output, result): + # A hierarchy's entry goes with it: on its metadata, which the copy shares. _set_output(self, output, result) def _tree(self, node): @@ -661,10 +894,15 @@ def _tree(self, node): return out def _unit(self, vtk_data): - """The filter's result for one unit: a dataset or a partitioned dataset.""" + """The filter's result for one unit: a dataset, a partitioned dataset or an AMR + hierarchy.""" + from tack.dataset import AMRDataSet from tack.interop.vtk import dataset_to_vtk - return dataset_to_vtk(self.execute(self._dataset(vtk_data, self.device)), copy=False) + result = self.execute(self._dataset(vtk_data, self.device)) + if isinstance(result, AMRDataSet): + return _amr_output(vtk_data, result) + return dataset_to_vtk(result, copy=False) # ── The filter ────────────────────────────────────────────────── @@ -713,6 +951,29 @@ class _FieldAlgorithm(TackAlgorithm): OUTPUT = None + def _amr_output_type(self): + return vtkOverlappingAMR + + def _by_box(self, vtk_amr): + """The hierarchy with the fields computed on each box, its blocks shallow copies + of the input's with the new arrays; its Tack hierarchy the input's with them, + on the same boxes and marks.""" + from tack.dataset import AMRDataSet + + amr = self._dataset(vtk_amr, self.device) + out = _amr_like(vtk_amr) + levels = [[] for _ in range(amr.num_levels)] + whole = True + for (level, index, block), (_, box) in zip(_blocks(vtk_amr), amr): + one, kept = self._with_fields(block, box) + out.SetDataSet(level, index, one) + levels[level].append(kept) + whole &= kept is not None + if whole: + _remember(out, AMRDataSet._grown(amr, levels, amr._extents, amr.layers, + amr._fills)) + return out + def _unit(self, vtk_data): data = self._dataset(vtk_data, self.device) if isinstance(vtk_data, vtkPartitionedDataSet): @@ -737,7 +998,7 @@ def _set_unit_output(self, output, result): output.CompositeShallowCopy(result) return output.ShallowCopy(result) - kept = _datasets.get(result.GetAddressAsString("vtkObject")) + kept = _entry(result) if kept is not None: # The output is a new object over the same arrays: the same Tack dataset. _remember(output, kept[1][None]) @@ -903,6 +1164,8 @@ class _Seeded(_Whole): INPUTS = 2 REQUIRED_INPUTS = 1 + # Box by box, a path would end where a finer box covers the coarse one. + BY_BOX = False def SetSourceConnection(self, port): self.SetInputConnection(1, port) @@ -960,12 +1223,14 @@ class Probe(TackAlgorithm): input port 1 (``SetSourceConnection``, ``SetSourceData``) -- at the points of the input (port 0), as point data on the input's structure, with ``valid`` (0 where no source cell holds the point). A partitioned source is probed as a whole, and a - tree as the partitioned dataset of its blocks; a tree on port 0 keeps its - structure. It runs on the source's device unless ``device`` says otherwise, the - input moved there.""" + tree as the partitioned dataset of its blocks; an AMR hierarchy at each point's + finest box (as ``vtkAMRResampleFilter``, onto any points); a tree on port 0 keeps + its structure. It runs on the source's device unless ``device`` says otherwise, + the input moved there.""" OUTPUT = "structure" INPUTS = 2 + AMR_PORTS = (1,) FUNCTION = "probe" EXCLUDE = ("where",) @@ -981,7 +1246,7 @@ def _prepare(self, inInfo): import tack.dataset as td source = self.GetInputData(inInfo, 1, 0) - if isinstance(source, _TREES): + if isinstance(source, _TREES) and not isinstance(source, vtkOverlappingAMR): # Every dataset in the tree, as the pieces of one partitioned dataset. leaves = vtkPartitionedDataSet() iterator = source.NewIterator() @@ -994,7 +1259,16 @@ def _prepare(self, inInfo): self._source = self._dataset(source, self.device) self._where_device = self.device if self._where_device is None and not isinstance(self._source, td.PartitionedDataSet): - self._where_device = self._source.device + self._where_device = _device_of(self._source) + + def _level(self, pieces): + import tack.dataset as td + + arguments = self.arguments() + return td.PartitionedDataSet( + [td.probe(self._source, piece if self._where_device is None + else piece.to_device(self._where_device), **arguments) + for piece in pieces], world=pieces.world) def _unit(self, vtk_data): import tack.dataset as td @@ -1512,9 +1786,10 @@ class UpwindFlux(_OnFaces): class Voxelize(_Whole): """``td.voxelize``: the input on the voxels of ``grid`` (a ``vtkImageData`` or a Tack image), with ``fraction``, the share of each voxel the input covers, and the - fields' means.""" + fields' means. An AMR hierarchy is sampled at each sample's finest box.""" OUTPUT = "vtkImageData" + AMR_PORTS = (0,) FUNCTION = "voxelize" DATASETS = ("grid",) @@ -1672,9 +1947,9 @@ class _OverPieces(TackAlgorithm): whole (a single dataset is one piece).""" def execute(self, data): - from tack.dataset import PartitionedDataSet + from tack.dataset import AMRDataSet, PartitionedDataSet - if not isinstance(data, PartitionedDataSet): + if not isinstance(data, (PartitionedDataSet, AMRDataSet)): # This process's one piece, beside the other processes' (``world``). data = PartitionedDataSet([data], world=self._run_world) return self.over_pieces(data) @@ -1683,13 +1958,18 @@ def over_pieces(self, pieces): return self._function()(pieces, **self.arguments()) def _output_type(self, input): - return vtkPartitionedDataSet + return vtkOverlappingAMR if isinstance(input, vtkOverlappingAMR) else vtkPartitionedDataSet class GhostCells(_OverPieces): """``td.ghost_cells``: each piece grown by ``layers`` of its neighbours' cells, - marked ``DUPLICATECELL``. VTK's ``vtkGhostCellsGenerator``.""" + marked ``DUPLICATECELL``. VTK's ``vtkGhostCellsGenerator``. An AMR hierarchy's + boxes each grow within their level, a ghost a sibling holds copied from it, else + from the coarser level (AMReX's ``FillBoundary`` and ``FillPatchTwoLevels``): a + hierarchy of the input's boxes, its blocks the grown ones, which the AMR filters + after it take as they are.""" + AMR_PORTS = (0,) FUNCTION = "ghost_cells" layers = Parameter(1, "the layers of ghost cells") @@ -1702,6 +1982,7 @@ class GenerateGlobalIds(_OverPieces): VTK's ``vtkGenerateGlobalIds``.""" FUNCTION = "generate_global_ids" + BY_BOX = False # a level's ids would repeat another's tolerance = Parameter(0.0, "how far apart two pieces' copies of a point may be") @@ -1744,16 +2025,66 @@ class GhostPoints(_OverPieces): FUNCTION = "ghost_points" +# ── AMR hierarchies ───────────────────────────────────────────────── + +class AverageDown(TackAlgorithm): + """``td.average_down``: an AMR hierarchy with its finer levels' values of ``field`` in + the coarse cells and points they cover, finest level first -- a refined cell the + volume-weighted mean of the fine cells over it, a refined point the fine point's + value at its location. AMReX's ``average_down``. The input stays as it was: the + output's ``field`` is a copy, its other arrays the input's.""" + + OUTPUT = None + AMR_PORTS = (0,) + FUNCTION = "average_down" + FIELDS = ("field",) + + def compute(self, data): + from tack.dataset import AMRDataSet + + if not isinstance(data, AMRDataSet): + raise TypeError("AverageDown takes an AMR hierarchy (vtkOverlappingAMR), not a " + f"{type(data).__name__}") + field = self.arguments()["field"] + out = _with_copy_of(data, field) + self._function()(out, field) # in place + return out + + +def _with_copy_of(amr, key): + """``amr`` with each box's field ``key`` a copy of it, the rest shared.""" + import tack + from tack.dataset import AMRDataSet, DataSet, Field + + levels = [] + for pieces in amr.levels: + boxes = [] + for box in pieces: + field = box.fields[key] + values = field.values + # An implicit array stays: average_down refuses it, naming it. + copied = Field(field.space, values.copy() if isinstance(values, tack.Array) + else values, field.unit) + fields = {name: copied if f is field else f for name, f in box.fields.items()} + boxes.append(DataSet(box.topology, box.geometry, fields=fields, sets=box.sets, + ghosts=box.ghosts, global_ids=box.global_ids, + **box._frame())) + levels.append(boxes) + return AMRDataSet._grown(amr, levels, amr._extents, amr.layers, amr._fills) + + # ── Reductions: one answer for all the input ──────────────────────── class _Reduction(TackAlgorithm): """A filter whose answer is one result for all its input: every dataset in it (a partitioned dataset's pieces, a tree's datasets) on every rank, reduced together, as VTK's ``vtkDescriptiveStatistics`` and ``vtkIntegrateAttributes`` reduce - composite data. Every rank gets the result. ``reduce(pieces)`` gives the VTK + composite data; an AMR hierarchy's boxes with their marks, so what a finer box + covers and what a sibling holds count once. Every rank gets the result. ``reduce(pieces)`` gives the VTK output from a ``PartitionedDataSet`` of every dataset in the input.""" OUTPUT = "vtkTable" + AMR_PORTS = (0,) fields = Parameter(None, "the fields: names, or (name, 'points' / 'cells'); by " "default every field. A vector field is reduced a " @@ -1773,7 +2104,12 @@ def RequestData(self, request, inInfo, outInfo): from tack.dataset import PartitionedDataSet self._run_world = self._world() + self._check_inputs(inInfo) input = self.GetInputData(inInfo, 0, 0) + if isinstance(input, vtkOverlappingAMR): + hierarchy = _as_dataset(input, self.device, self._run_world) + self.GetOutputData(outInfo, 0).ShallowCopy(self.reduce(hierarchy._as_pieces())) + return 1 if isinstance(input, (*_TREES, vtkPartitionedDataSet)): leaves = [] iterator = input.NewIterator() @@ -1999,4 +2335,494 @@ def reduce(self, pieces): return grid + +# ── Temporal filters ──────────────────────────────────────────────── + +class _Temporal(TackAlgorithm): + """A temporal filter of ``tack.dataset.temporal`` -- an object fed one step at a + time -- in a VTK pipeline, in one of two modes (``mode``): + + - **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. An + execution at the time of the last step gives that output again, without + taking the step twice; a step back in time, a parameter changed or VTK's + ``NO_PRIOR_TEMPORAL_ACCESS`` saying reset start over. + - **Driving** the filter asks its input for the steps it needs, one per pass + (``CONTINUE_EXECUTING``), never one it has taken since it last started over. + A filter that accumulates (``ACCUMULATES``: statistics, particles, paths) + takes every step (or those in ``time_range``) and gives one output after the + last, with no time: time is in its data, and it reports no time steps + downstream. One whose answer is at a time (streaklines, a derivative) gives + it at the time requested of it, going on from the steps it has taken; a + request before them starts over. + + By default a filter drives where its input has time steps and VTK does not say + the pipeline cannot go back in time (``NO_PRIOR_TEMPORAL_ACCESS``, which an in + situ adaptor sets). Driving to a step gives the output the same steps give in + situ. A tree or an AMR hierarchy is refused; a partitioned dataset's pieces each + keep their state, as Tack's filters do. Subclasses give ``_make()``, the Tack + filter, and ``_output(result, input, data, final)``, the VTK output for its + result.""" + + INPUTS = 1 + BY_BOX = False + #: Whether the filter accumulates over the steps (one output, after the last, + #: when driving), rather than answering at a time. + ACCUMULATES = True + #: Whether an answer at step k needs the steps before it (streaklines), or only + #: step k - 1 (a derivative). + NEEDS_EVERY_STEP = True + + mode = Parameter(None, "'in situ' or 'driving'; by default driving where the input " + "has time steps and the pipeline can go back in time") + + def __init__(self, **parameters): + super().__init__(**parameters) + self._driving = False + self._steps = [] + self._queue = [] + self._looping = False + self._start_over() + + def _start_over(self): + self._filter = None # the Tack filter, made at the first step + self._fed = None # driving: the index of the last step taken + self._time = None # in situ: the time of the last step taken + self._made_at = None # the parameter changes it was made at + self._kept = None # (output, index or time) of the last output + self._results = [] # driving, accumulating: what each step gave + self._stale = False + + def _changes(self): + return self.__dict__.get("_parameter_changes", 0) + + def _drives(self, info): + """Whether this run drives the pipeline (``mode``, else as the input says).""" + s = _pipeline() + mode = (self.mode or "").replace("_", " ") + if mode in ("in situ", "driving"): + return mode == "driving" + if info.Has(s.NO_PRIOR_TEMPORAL_ACCESS()): + return False + return info.Has(s.TIME_STEPS()) and info.Length(s.TIME_STEPS()) > 0 + + # ── The pipeline ──────────────────────────────────────────────── + + def RequestInformation(self, request, inInfo, outInfo): + s = _pipeline() + info = inInfo[0].GetInformationObject(0) + self._driving = self._drives(info) + self._steps = list(info.Get(s.TIME_STEPS())) if info.Has(s.TIME_STEPS()) else [] + # The pipeline before this filter, or the filter, changed: what was taken + # of the steps is stale (in situ, every new step is such a change, and + # the steps' times say what is new instead). + self._stale = True + if self._driving and self.ACCUMULATES: + out = outInfo.GetInformationObject(0) + out.Remove(s.TIME_STEPS()) + out.Remove(s.TIME_RANGE()) + return 1 + + def RequestUpdateExtent(self, request, inInfo, outInfo): + if not self._driving: + return 1 # the step the input holds + s = _pipeline() + if not self._looping: + self._plan(outInfo.GetInformationObject(0)) + info = inInfo[0].GetInformationObject(0) + if self._queue: + info.Set(s.UPDATE_TIME_STEP(), self._steps[self._queue[0]]) + elif self._fed is not None: + # Nothing to take: the input keeps the step it gave last, unexecuted. + info.Set(s.UPDATE_TIME_STEP(), self._steps[self._fed]) + return 1 + + def _plan(self, out): + """The steps this run takes, in order (``_queue``): none where the output kept + is the answer.""" + s = _pipeline() + stale = self._stale or self._made_at != self._changes() + self._stale = False + if self.ACCUMULATES: + window = self.time_range + steps = [i for i, t in enumerate(self._steps) + if window is None or window[0] <= t <= window[1]] + if not steps: + raise ValueError(f"{type(self).__name__}: no time step in {window}") + if stale or self._kept is None: + self._start_over() + self._queue = steps + else: + self._queue = [] + return + wanted = out.Get(s.UPDATE_TIME_STEP()) if out.Has(s.UPDATE_TIME_STEP()) else None + k = 0 if wanted is None else max( + [i for i, t in enumerate(self._steps) if t <= wanted] or [0]) + fed = None if stale else self._fed + if fed is not None and self._kept is not None and self._kept[1] == k: + self._queue = [] # the answer is kept + elif fed is not None and (fed < k if self.NEEDS_EVERY_STEP else fed == k - 1): + self._queue = list(range(fed + 1, k + 1)) # on from the steps taken + else: + self._start_over() + first = 0 if self.NEEDS_EVERY_STEP else max(k - 1, 0) + self._queue = list(range(first, k + 1)) + + def RequestData(self, request, inInfo, outInfo): + s = _pipeline() + self._run_world = self._world() + self._check_inputs(inInfo) + input = self.GetInputData(inInfo, 0, 0) + if isinstance(input, _TREES): + raise TypeError(f"{type(self).__name__} takes a dataset or a partitioned " + f"dataset, not a {input.GetClassName()}") + output = self.GetOutputData(outInfo, 0) + try: + if self._driving: + self._drive(request, input, output) + else: + self._in_situ(inInfo[0].GetInformationObject(0), input, output) + except BaseException: + request.Remove(s.CONTINUE_EXECUTING()) + self._looping = False + self._start_over() + raise + return 1 + + def _drive(self, request, input, output): + s = _pipeline() + if not self._queue: + self._set(output, *self._kept) + return + k = self._queue.pop(0) + data = self._dataset(input, self.device) + if self._filter is None: + self._filter = self._make(data) + self._made_at = self._changes() + result = self._filter.update(data) + self._fed = k + if self.ACCUMULATES: + self._results.append((k, result)) + if self._queue: + request.Set(s.CONTINUE_EXECUTING(), 1) + self._looping = True + return + request.Remove(s.CONTINUE_EXECUTING()) + self._looping = False + if self.ACCUMULATES: + self._kept = (self._output(self._results, input, data, final=True), None) + self._results = [] + else: + self._kept = (self._output(result, input, data, final=False), k) + self._set(output, *self._kept) + + def _in_situ(self, info, input, output): + from tack.interop.vtk import _time_of + + s = _pipeline() + data = self._dataset(input, self.device) + t = _time_of(data) + if t is None: + raise ValueError(f"{type(self).__name__} needs each step's time " + "(DATA_TIME_STEP)") + reset = (self._filter is None or self._made_at != self._changes() + or (self._time is not None and t < self._time) + or (info.Has(s.NO_PRIOR_TEMPORAL_ACCESS()) + and info.Get(s.NO_PRIOR_TEMPORAL_ACCESS()) + == s.NO_PRIOR_TEMPORAL_ACCESS_RESET)) + if not reset and t == self._time: + self._set(output, *self._kept) # no new step + return + if reset: + self._start_over() + self._filter = self._make(data) + self._made_at = self._changes() + result = self._filter.update(data) + self._time = t + self._kept = (self._output(result, input, data, final=False), t) + self._set(output, *self._kept) + + def _set(self, output, result, at): + """Make ``output`` the VTK ``result``: at its step's time, or at none.""" + _set_output(self, output, result) + info = output.GetInformation() + if at is None: + info.Remove(vtkDataObject.DATA_TIME_STEP()) + else: + info.Set(vtkDataObject.DATA_TIME_STEP(), + self._steps[at] if self._driving else at) + + # ── What a subclass gives ─────────────────────────────────────── + + def _make(self, data): + raise NotImplementedError + + def _output(self, result, input, data, final): + """The VTK output: for ``result``, the Tack filter's for the step ``data`` (the + VTK ``input``); or when driving a filter that accumulates (``final``), for + ``result``, every step's ``(index, result)``.""" + raise NotImplementedError + + def _fields(self, data): + """``fields`` as Tack's keys; by default every field of the step.""" + from tack.dataset import PartitionedDataSet + + first = next(iter(data)) if isinstance(data, PartitionedDataSet) else data + if self.fields is None: + return list(first.fields) + chosen = [self.fields] if isinstance(self.fields, (str, tuple)) else self.fields + return [_field_argument(value) for value in chosen] + + +def _pipeline(): + from vtkmodules.vtkCommonExecutionModel import vtkStreamingDemandDrivenPipeline + + return vtkStreamingDemandDrivenPipeline + + +def _arrays_onto(vtk_data, pieces, names, keep): + """A shallow copy of the VTK dataset or partitioned dataset ``vtk_data`` with the + fields ``names`` of the matching Tack ``pieces`` (a dataset, or a partitioned one) + as arrays; its own arrays kept if ``keep``, else only its ghost marks.""" + from tack.interop.vtk import _label, _set_values, _to_vtk_values + + def onto(block, piece): + out = block.NewInstance() + out.ShallowCopy(block) + if not keep: + for attributes in (out.GetPointData(), out.GetCellData()): + ghosts = attributes.GetArray("vtkGhostType") + attributes.Initialize() + if ghosts is not None: + attributes.AddArray(ghosts) + for name in names: + field = piece.fields[name] + attributes = out.point_data if _on_points(field.space) else out.cell_data + _set_values(attributes, name, _to_vtk_values(field.values, share=True, + name=name)) + _label(attributes[name], field.unit) + return out + + if isinstance(vtk_data, vtkPartitionedDataSet): + out = vtk_data.NewInstance() + for block, piece in zip(vtk_data, pieces): + out.append(onto(block, piece)) + return out + return onto(vtk_data, pieces) + + +def _count(vtk_data, name, value): + """``vtk_data`` (a dataset, or each partition) with the field data array ``name``, + one integer.""" + from vtkmodules.vtkCommonCore import vtkTypeInt64Array + + blocks = list(vtk_data) if isinstance(vtk_data, vtkPartitionedDataSet) else [vtk_data] + for block in blocks: + array = vtkTypeInt64Array() + array.SetName(name) + array.InsertNextValue(int(value)) + block.GetFieldData().AddArray(array) + return vtk_data + + +class RunningStatistics(_Temporal): + """``td.RunningStatistics``: per point or cell, ``_mean``, ``_variance`` + (the population variance), ``_min`` and ``_max`` of ``fields`` (every + field by default) over the steps, on the input's structure, with field data + ``count``, the steps taken. In situ, the statistics up to the step, at its time; + driving, over every step (or those in ``time_range``), at no time. The steps must + be on one mesh. VTK's ``vtkTemporalStatistics``.""" + + OUTPUT = None + + fields = Parameter(None, "the fields: names, or (name, 'points' / 'cells'); by " + "default every field") + time_range = Parameter(None, "driving: (first, last), the times of the steps to " + "take; by default every step") + + def _make(self, data): + import tack.dataset as td + + self._names = self._fields(data) + return td.RunningStatistics(self._names) + + def _output(self, result, input, data, final): + from tack.dataset.dataset import field_name + + if final: + count, result = len(result), result[-1][1] + else: + count = self._filter.count + names = [f"{field_name(name)}_{suffix}" for name in self._names + for suffix in ("mean", "variance", "min", "max")] + return _count(_arrays_onto(input, result, names, keep=False), "count", count) + + +class TimeDerivative(_Temporal): + """``td.TimeDerivative``: the step with ``_dt`` beside its arrays, each of + ``fields``' (every field by default) rate of change since the step before, a + backward difference. The first step has no step before: it goes out without + them. Driving, at the time requested, from that step and the one before it + (only the requested one where the one before is kept). VTK has none.""" + + OUTPUT = None + ACCUMULATES = False + NEEDS_EVERY_STEP = False + + fields = Parameter(None, "the fields: names, or (name, 'points' / 'cells'); by " + "default every field") + copy = Parameter(True, "keep a copy of the step before (False: references, where " + "the step's arrays are not overwritten by the next)") + + def _make(self, data): + import tack.dataset as td + + self._names = self._fields(data) + return td.TimeDerivative(self._names, copy=self.copy) + + def _output(self, result, input, data, final): + from tack.dataset.dataset import field_name + + if result is None: # the first step + return _arrays_onto(input, data, [], keep=True) + return _arrays_onto(input, result, [f"{field_name(name)}_dt" for name in self._names], + keep=True) + + +class _Particles(_Temporal): + """A temporal filter of particles from seeds: ``seeds`` (an ``(n, 3)`` array of + world positions) or the points of the data on port 1, read when it starts.""" + + INPUTS = 2 + REQUIRED_INPUTS = 1 + OUTPUT = "vtkUnstructuredGrid" + + field = Parameter(None, "the velocity field (required): a name, or (name, 'points' " + "/ 'cells')") + seeds = Parameter(None, "the seeds, an (n, 3) array of world positions; by default " + "the points of the data on port 1") + integrator = Parameter("rk4", "rk4 or euler") + step_size = Parameter(None, "the longest substep, in the steps' time unit; by default " + "one substep from step to step") + copy = Parameter(True, "keep a copy of the step before (False: references)") + + def SetSourceConnection(self, port): + self.SetInputConnection(1, port) + + def SetSourceData(self, data): + self.SetInputDataObject(1, data) + + def RequestData(self, request, inInfo, outInfo): + self._source = (self.GetInputData(inInfo, 1, 0) + if inInfo[1].GetNumberOfInformationObjects() else None) + return super().RequestData(request, inInfo, outInfo) + + def _arguments(self): + field = self.field if self.field is not None else self._input_array(0) + if field is None: + raise ValueError(f"{type(self).__name__} needs 'field', the velocity") + seeds = self.seeds + if seeds is None and self._source is not None: + seeds = _points_of(self._source) + if seeds is None: + raise ValueError(f"{type(self).__name__} needs seeds: 'seeds', or data on " + "port 1") + return {"field": _field_argument(field), "seeds": seeds, + "integrator": self.integrator, "step_size": self.step_size, + "copy": self.copy} + + def _output(self, result, input, data, final): + from tack.interop.vtk import dataset_to_vtk + + out = dataset_to_vtk(result, copy=False) + # Of pieces, this rank's paths or lines are one dataset: its one partition. + return _as_partitioned(out) if isinstance(input, vtkPartitionedDataSet) else out + + +class Advection(_Particles): + """``td.Advection``: particles from ``seeds`` carried by the velocity ``field`` + through the steps, the velocity linear in time between two steps. In situ, the + particles at the step's time: a point cloud with ``seed``, ``source``, ``release``, + ``time`` and ``status``. Driving, the particles at every step taken (or those in + ``time_range``) in one point cloud, with ``step`` (its index among the input's + time steps) beside those, at no time. Over MPI ranks each rank releases the seeds + it is given.""" + + time_range = Parameter(None, "driving: (first, last), the times of the steps to " + "take; by default every step") + + def _make(self, data): + import tack.dataset as td + + return td.Advection(**self._arguments()) + + def _output(self, result, input, data, final): + if not final: + return super()._output(result, input, data, final) + return super()._output(_every_step(result), input, data, final) + + +def _every_step(results): + """One point cloud of every step's particles (``(index, cloud)``, each a cloud or a + partitioned dataset of this rank's), with ``step``, at no time.""" + import numpy as np + + import tack.dataset as td + + clouds = [] + for k, cloud in results: + for piece in (cloud if isinstance(cloud, td.PartitionedDataSet) else [cloud]): + clouds.append((k, piece)) + first = clouds[0][1] + names = ("seed", "source", "release", "time", "status") + fields = {name: np.concatenate([np.asarray(c.fields[name].values.numpy()) + for _, c in clouds]) for name in names} + fields["step"] = np.concatenate([np.full(c.num_points, k, np.int64) for k, c in clouds]) + unit = first.fields["time"].unit + out = td.point_cloud(np.concatenate([c.world_positions() for _, c in clouds]), + fields=fields, origin=first.origin, device=first.device, + unit=first.unit, + units={"release": unit, "time": unit} if unit else None) + last = results[-1][1] + if isinstance(last, td.PartitionedDataSet): + return td.PartitionedDataSet([out], world=last.world) + return out + + +class Pathlines(_Particles): + """``td.Pathlines``: each particle's path from ``seeds`` through the velocity + ``field``: a polyline per particle, ``seed`` and ``time`` on its points. In situ, + the paths up to the step, at its time; driving, through every step taken (or + those in ``time_range``), at no time. VTK's ``vtkParticlePathFilter``.""" + + time_range = Parameter(None, "driving: (first, last), the times of the steps to " + "take; by default every step") + + def _make(self, data): + import tack.dataset as td + + return td.Pathlines(**self._arguments()) + + def _output(self, result, input, data, final): + return super()._output(result[-1][1] if final else result, input, data, final) + + +class Streaklines(_Particles): + """``td.Streaklines``: the streaklines at a time of particles released at ``seeds`` + through the velocity ``field``, ``releases`` at each seed from step to step: a + polyline per seed through its particles still in the mesh, newest first, with + ``source``, ``release`` and ``time`` on the points, at the step's time. Driving, + at the time requested, going on from the steps taken. VTK's + ``vtkStreaklineFilter``.""" + + ACCUMULATES = False + + releases = Parameter(1, "the particles released at each seed from step to step") + + def _make(self, data): + import tack.dataset as td + + return td.Streaklines(**self._arguments(), releases=self.releases) + + del tack diff --git a/packages/tack-vis/tests/_temporal_mpi.py b/packages/tack-vis/tests/_temporal_mpi.py index d1cd3c3c5..7d820cd5b 100644 --- a/packages/tack-vis/tests/_temporal_mpi.py +++ b/packages/tack-vis/tests/_temporal_mpi.py @@ -18,10 +18,11 @@ def main(): comm = MPI.COMM_WORLD - got, want = case.advect_across(tc.World(comm)) - np.testing.assert_allclose(got.world_positions(), want.world_positions(), atol=1e-5) - np.testing.assert_array_equal(got.fields["status"].values.numpy(), - want.fields["status"].values.numpy()) + for uneven in (False, True): + got, want = case.advect_across(tc.World(comm), uneven) + np.testing.assert_allclose(got.world_positions(), want.world_positions(), atol=1e-5) + np.testing.assert_array_equal(got.fields["status"].values.numpy(), + want.fields["status"].values.numpy()) done = comm.gather(True, root=0) if comm.rank == 0 and all(done): print(f"ok: {comm.size} ranks", flush=True) diff --git a/packages/tack-vis/tests/_vtk_algorithms_mpi.py b/packages/tack-vis/tests/_vtk_algorithms_mpi.py index fd73f6fbb..a21434cdf 100644 --- a/packages/tack-vis/tests/_vtk_algorithms_mpi.py +++ b/packages/tack-vis/tests/_vtk_algorithms_mpi.py @@ -1,10 +1,15 @@ +# SPDX-FileCopyrightText: Copyright (c) Kitware Inc. +# SPDX-License-Identifier: BSD-3-Clause + """``test_vtk_algorithms``' check across MPI ranks: ``mpiexec -n 3 python _vtk_algorithms_mpi.py comm|vtk``. Each rank holds one slab of VTK's wavelet, the slabs sharing their planes of points, in a ``vtkPartitionedDataSet`` of one -partition; the filters over pieces run in a VTK pipeline on every rank. The world -is given as ``world=`` (an mpi4py communicator: ``comm``) or taken from VTK's global -controller (``vtk``, which needs a VTK built with MPI and vtkParallelMPI4Py). -Rank 0 prints ``ok``, or the failure.""" +partition; the filters over pieces run in a VTK pipeline on every rank. Then each +holds blocks of VTK's AMR Gaussian pulse, the others' empty, as VTK spreads a +hierarchy; the AMR filters give the whole's. The world is given as ``world=`` (an +mpi4py communicator: ``comm``) or taken from VTK's global controller (``vtk``, +which needs a VTK built with MPI and vtkParallelMPI4Py). Rank 0 prints ``ok``, or +the failure.""" import sys @@ -72,8 +77,109 @@ def main(): merged = ta.MergePieces(**options)(pieces) if rank == 0: assert merged.number_of_points == points, (merged.number_of_points, points) + hierarchy(comm, options) + temporal(comm, options) + if rank == 0: print(f"ok: {size} ranks, world from {mode}", flush=True) +def flow(t): + """A step at ``t`` as VTK's: an image with a velocity ``v`` along x, turning.""" + import tack.dataset as td + from tack.interop.vtk import dataset_to_vtk + + grid = td.uniform_grid((13, 7, 3)) + grid.time = t + v = np.zeros((grid.num_points, 3), np.float32) + v[:, 0], v[:, 1] = 1.0, 0.2 * t + grid.fields["v"] = td.Field(td.H1(grid), tack.array(v, tack.vec3())) + return dataset_to_vtk(grid) + + +def temporal(comm, options): + """Pathlines in situ over slabs of each step, the seeds on rank 0: particles cross + from rank to rank, and the paths are the whole's.""" + from tack.communication import World + + rank, size = comm.Get_rank(), comm.Get_size() + seeds = np.array([[0.5, 2.0, 1.0], [1.5, 4.0, 1.0]]) + cuts = np.linspace(0, 12, size + 1).round().astype(int) + spread = ta.Pathlines(field="v", seeds=seeds if rank == 0 else np.zeros((0, 3)), + **options) + whole = ta.Pathlines(field="v", seeds=seeds, world=World()) + for t in (0.0, 1.0, 2.0, 3.0): + step = flow(t) + slab = dm.vtkPartitionedDataSet() + slab.append(vtkExtractVOI(voi=(int(cuts[rank]), int(cuts[rank + 1]), 0, 6, 0, 2))(step)) + slab.GetInformation().Set(dm.vtkDataObject.DATA_TIME_STEP(), t) + spread.SetInputDataObject(0, slab) + spread.Update() + assert spread.error is None, spread.error + whole.SetInputDataObject(0, step) + whole.Update() + mine = spread.GetOutputDataObject(0) + points = [np.asarray(block.points.data) for block in mine if block.number_of_points] + every = comm.gather(np.concatenate(points) if points else np.zeros((0, 3))) + if rank == 0: + got = np.concatenate(every) + expected = np.asarray(whole.GetOutputDataObject(0).points.data) + assert got.shape == expected.shape, (got.shape, expected.shape) + np.testing.assert_allclose(got, expected, rtol=1e-6, atol=1e-6) + + +def pulse(): + from vtkmodules.vtkFiltersAMR import vtkAMRGaussianPulseSource + + return vtkAMRGaussianPulseSource(dimension=3)() + + +def hierarchy(comm, options): + """The pulse's blocks dealt to the ranks in turn: the statistics and integral of the + spread hierarchy, its averages down and its ghost cells, and a filter box by box, + are the whole's.""" + from tack.communication import World + + rank, size = comm.Get_rank(), comm.Get_size() + + alone = {"world": World()} + whole, spread = pulse(), pulse() + mine = [] + blocks = [(level, index) for level in range(spread.GetNumberOfLevels()) + for index in range(spread.GetNumberOfBlocks(level))] + for k, (level, index) in enumerate(blocks): + if k % size == rank: + mine.append((level, index)) + else: + # Through the level's partitioned dataset: SetDataSet ignores None. + spread.GetPartitionedDataSet(level).SetPartition(index, None) + for make in (lambda **o: ta.Statistics(fields=["Gaussian-Pulse"], **o), + lambda **o: ta.GhostCells(layers=1, **o) >> ta.Statistics( + fields=["Gaussian-Pulse"], **o)): + expected, got = make(**alone)(whole), make(**options)(spread) + for name in ("n", "sum", "min", "max"): + a = got.GetColumnByName(name).GetValue(0) + b = expected.GetColumnByName(name).GetValue(0) + assert abs(a - b) <= 1e-9 * max(1.0, abs(b)), (name, a, b) + def cells(collection, level): + return sum(block.number_of_cells for block in collection.GetPartitionedDataSet(level)) + + threshold = {"field": ("Gaussian-Pulse", "cells"), "lower": 0.0, "upper": 1.0} + kept = ta.Threshold(**threshold, **options)(spread) # box by box + whole_kept = ta.Threshold(**threshold, **alone)(whole) + for level in range(whole.GetNumberOfLevels()): + total = comm.allreduce(cells(kept, level)) + assert total == cells(whole_kept, level), (level, total) + volume = ta.Integrate(**options)(spread).cell_data["Volume"].GetValue(0) + assert abs(volume - ta.Integrate(**alone)(whole).cell_data["Volume"].GetValue(0)) < 1e-12 + averaged = ta.AverageDown(field="Gaussian-Pulse", **options)(spread) + reference = ta.AverageDown(field="Gaussian-Pulse", **alone)(whole) + for level, index in mine: + got = np.asarray(averaged.GetDataSetAsCartesianGrid(level, index) + .cell_data["Gaussian-Pulse"]) + want = np.asarray(reference.GetDataSetAsCartesianGrid(level, index) + .cell_data["Gaussian-Pulse"]) + assert np.array_equal(got, want), (level, index) + + if __name__ == "__main__": main() diff --git a/packages/tack-vis/tests/test_temporal.py b/packages/tack-vis/tests/test_temporal.py index e8d2d2808..fff56b5cb 100644 --- a/packages/tack-vis/tests/test_temporal.py +++ b/packages/tack-vis/tests/test_temporal.py @@ -512,13 +512,24 @@ def _swirl(x, t): return v -def advect_across(world): +def advect_across(world, uneven=False): """This rank's particles (its own seeds) advected through steps whose pieces are - shared out over ``world``'s ranks, and the same seeds through the whole steps.""" + shared out over ``world``'s ranks, and the same seeds through the whole steps. + ``uneven``: every seed on rank 0, one leaving the mesh, and a release on the + last rank alone, so ranks have no particles, stop at other substeps and + release at other steps, while the probes are collective.""" seeds = np.array([[7.0, 5.0, 1.0], [5.0, 3.5, 1.0], [4.0, 6.0, 1.0], - [6.0, 6.5, 1.0]])[world.rank::world.size] + [6.0, 6.5, 1.0]]) + if uneven: + seeds = np.concatenate([seeds, [[9.6, 5.0, 1.0]]]) if world.rank == 0 else \ + np.zeros((0, 3)) + else: + seeds = seeds[world.rank::world.size] shared, whole = td.Advection("v", seeds, step_size=0.1), td.Advection("v", seeds, step_size=0.1) + if uneven and world.rank == world.size - 1: + for advection in (shared, whole): + advection.inject([[5.0, 5.0, 1.0]], time=0.25) for t in (0.0, 0.5, 1.0): step = _flow(t, _swirl) pieces = list(td.split(step, (2, 2, 1)))[world.rank::world.size] diff --git a/packages/tack-vis/tests/test_vtk_algorithms.py b/packages/tack-vis/tests/test_vtk_algorithms.py index edaecddf1..1ce5aa57d 100644 --- a/packages/tack-vis/tests/test_vtk_algorithms.py +++ b/packages/tack-vis/tests/test_vtk_algorithms.py @@ -242,18 +242,35 @@ def test_a_changed_dataset_is_converted_again(): def test_an_entry_goes_with_its_vtk_dataset(): + """The entry is an attribute of the dataset's wrapper: it outlives the wrapper while + the dataset lives (a fresh wrapper has it) and goes with the dataset.""" import gc + import weakref - grid = _rotation() - ta._as_dataset(grid) - address = grid.GetAddressAsString("vtkObject") - assert address in ta._datasets - del grid + holder = vtk.vtkPartitionedDataSet() + holder.append(_rotation()) + data = weakref.ref(ta._as_dataset(holder.GetPartition(0))) + gc.collect() + assert data() is not None and ta._as_dataset(holder.GetPartition(0)) is data() + holder.SetPartition(0, None) gc.collect() - # The deletion is noted, and the entry goes at the cache's next use. - assert address in ta._deleted - ta._drop_deleted() - assert address not in ta._datasets and address not in ta._observers + assert data() is None + + +def test_an_entry_is_stale_after_modified_or_a_new_time(): + """``Modified()`` (here an array's) or another time step converts again; values + written in place without ``Modified()`` reach the fields, which share VTK's memory, + but the entry stays, as VTK's own pipeline would not run again either.""" + grid = _rotation() + data = ta._as_dataset(grid) + np.asarray(grid.GetPointData().GetArray("v"))[0] = 7.0 + assert ta._as_dataset(grid) is data + assert data.fields["v"].values.numpy()[0, 0] == 7.0 + grid.GetPointData().GetArray("v").Modified() + moved = ta._as_dataset(grid) + assert moved is not data + grid.GetInformation().Set(vtk.vtkDataObject.DATA_TIME_STEP(), 3.0) + assert ta._as_dataset(grid) is not moved def test_a_field_filters_output_keeps_its_inputs_topology(): @@ -936,17 +953,16 @@ def test_each_wrapper_gives_its_tack_filters_output(name): _assert_same(_contents(out), _contents(dataset_to_vtk(direct(_tack(data()))))) -def test_every_algorithm_but_the_tables_and_amr_has_a_wrapper(): - """Every algorithm of the reference page, but those for later (statistics, - entropy, integrate, average_down), the in-place update_ghosts, and the finite - volume operators that take values on faces, which VTK keeps nowhere.""" +def test_every_algorithm_has_a_wrapper(): + """Every algorithm of the reference page, but the in-place update_ghosts and the + finite volume operators that take values on faces, which VTK keeps nowhere (and + the first wrappers, named otherwise).""" import re from pathlib import Path page = (Path(__file__).parents[3] / "docs/vis-guide/reference/algorithms.md").read_text() names = set(re.findall(r"^\| \[`(\w+)`\]", page, re.M)) - later = {"statistics", "entropy", "integrate", "average_down", "update_ghosts", - "net_flux", "divergence", "perot", "contour", "slice_plane", "clip", "threshold", + later = {"update_ghosts", "net_flux", "divergence", "perot", "contour", "slice_plane", "clip", "threshold", "extract_surface", "tetrahedralize", "surface_normals", "connected_components", "streamlines", "probe", "gradients", "point_elevation", "cell_quality", "ghost_cells", "generate_global_ids", "merge_pieces"} @@ -1047,28 +1063,28 @@ def ends(data): np.testing.assert_allclose(furthest_halves, furthest_whole, rtol=1e-6) -def test_a_dataset_the_collector_frees_leaves_its_observer_callable(): +def test_a_dataset_the_collector_frees_takes_its_entry(): """A dataset an algorithm read, freed by Python's collector in a cycle made after - it was read (a traceback's frames hold one so): its DeleteEvent observer, referenced - only by it, was cleared as garbage first, and calling it crashed Python (VTK bugs - 19). The observers are kept reachable, so this collects without crashing.""" + it was read (as a traceback's frames hold one): a DeleteEvent observer, referenced + only by it, was cleared as garbage first and calling it crashed Python (VTK bugs + 19). The entry is the wrapper's attribute, no observer: it goes with the dataset.""" import gc + import weakref class Holder: pass - addresses = [] + entries = [] for _ in range(3): data = _wavelet((-4, 4, -4, 4, -4, 4))() assert ta.Contour(field="RTData", value=150.0)(data).number_of_cells > 0 - addresses.append(data.GetAddressAsString("vtkObject")) + entries.append(weakref.ref(ta._as_dataset(data))) holder = Holder() holder.cycle = holder holder.data = data del data, holder gc.collect() - ta._drop_deleted() - assert not any(address in ta._observers for address in addresses) + assert all(entry() is None for entry in entries) # ── Reductions: Statistics, Entropy, Integrate ────────────────────── @@ -1161,3 +1177,504 @@ def test_a_surfaces_measure_is_its_area(): out = ta.Integrate()(surface) assert "Area" in out.cell_data.keys() assert out.cell_data["Area"].GetValue(0) == pytest.approx(6 * 16 * 16, rel=1e-12) + + +# ── AMR hierarchies ───────────────────────────────────────────────── + +def _pulse(): + """VTK's AMR Gaussian pulse: a level-0 box and two finer boxes in it, the pulse as + cell data.""" + from vtkmodules.vtkFiltersAMR import vtkAMRGaussianPulseSource + + return vtkAMRGaussianPulseSource(dimension=3)() + + +def _hierarchy(amr): + return vtk_to_dataset(amr, dtype=None) + + +def _amr_bounds(amr): + """``(lo, hi)``: a hierarchy's ``GetBounds()`` returns a pointer's name (9.7.1).""" + bounds = [0.0] * 6 + amr.GetBounds(bounds) + return np.asarray(bounds[::2]), np.asarray(bounds[1::2]) + + +PULSE = ("Gaussian-Pulse", td.Constant) + + +def test_an_amr_hierarchy_is_reduced_once_per_region(): + """What a finer box covers counts once: the statistics, entropy and integral of the + hierarchy are Tack's, the cells counted those no finer box covers.""" + amr = _pulse() + hierarchy = _hierarchy(amr) + stats = td.statistics(hierarchy, PULSE) + table = ta.Statistics(fields=["Gaussian-Pulse"])(amr) + for name in ta.Statistics.NAMES: + assert table.GetColumnByName(name).GetValue(0) == pytest.approx( + getattr(stats, name), rel=1e-12), name + every_cell = sum(amr.GetDataSetAsCartesianGrid(level, i).number_of_cells + for level in range(amr.GetNumberOfLevels()) + for i in range(amr.GetNumberOfBlocks(level))) + assert stats.n < every_cell + entropy = ta.Entropy(fields=["Gaussian-Pulse"], bins=16)(amr) + assert entropy.GetColumnByName("entropy").GetValue(0) == pytest.approx( + td.entropy(hierarchy, PULSE, bins=16), rel=1e-12) + integral = ta.Integrate(fields=["Gaussian-Pulse"])(amr) + lo, hi = _amr_bounds(amr) + assert integral.cell_data["Volume"].GetValue(0) == pytest.approx(np.prod(hi - lo), + rel=1e-12) + expected = td.integrate(hierarchy._as_pieces(), fields=[PULSE]).values + assert integral.cell_data["Gaussian-Pulse"].GetValue(0) == pytest.approx( + next(iter(expected.values())), rel=1e-12) + + +def test_probe_and_voxelize_take_an_amr_hierarchy_at_its_finest_boxes(): + amr = _pulse() + hierarchy = _hierarchy(amr) + lo, hi = _amr_bounds(amr) + image = vtk.vtkImageData(dimensions=(9, 9, 9), origin=tuple(lo), + spacing=tuple((hi - lo) / 8)) + probe = ta.Probe() + probe.SetInputDataObject(0, image) + probe.SetSourceData(amr) + probe.Update() + assert probe.error is None, probe.error + expected = td.probe(hierarchy, vtk_to_dataset(image, dtype=None)) + np.testing.assert_array_equal( + np.asarray(probe.GetOutputDataObject(0).point_data["Gaussian-Pulse"]), + expected.fields["Gaussian-Pulse"].values.numpy()) + voxels = ta.Voxelize(grid=image, fields=["Gaussian-Pulse"], samples=2)(amr) + assert isinstance(voxels, vtk.vtkImageData) + direct = td.voxelize(hierarchy, vtk_to_dataset(image, dtype=None), + fields=[PULSE], samples=2) + np.testing.assert_array_equal(np.asarray(voxels.cell_data["fraction"]), + direct.fields["fraction"].values.numpy()) + + +def test_ghost_cells_of_an_amr_hierarchy_feed_the_next_amr_filter(): + """The grown boxes go out as the input's hierarchy with larger blocks, and the next + filter takes the grown hierarchy as it was made, in a pipeline or given what + calling GhostCells returned (a shallow copy): VTK's blocks with ghost layers + could not be converted again. Its ghosts count nowhere.""" + amr = _pulse() + grown = ta.GhostCells(layers=1)(amr) + assert isinstance(grown, vtk.vtkOverlappingAMR) + expected = td.ghost_cells(_hierarchy(amr), 1) + for (level, box), (lv, i, block) in zip(expected, ta._blocks(grown)): + assert level == lv + assert block.GetDimensions() == tuple(box.topology.point_dims) + assert grown.GetAMRBox(lv, i).GetNumberOfCells() == amr.GetAMRBox( + lv, i).GetNumberOfCells() + assert ta._as_dataset(grown).layers == 1 + once = ta.Statistics(fields=["Gaussian-Pulse"])(amr) + pipeline = ta.GhostCells(layers=1) >> ta.Statistics(fields=["Gaussian-Pulse"]) + for table in (pipeline(amr), ta.Statistics(fields=["Gaussian-Pulse"])(grown)): + assert table.GetColumnByName("n").GetValue(0) == once.GetColumnByName( + "n").GetValue(0) + assert table.GetColumnByName("mean").GetValue(0) == pytest.approx( + once.GetColumnByName("mean").GetValue(0), rel=1e-12) + # The input's entry is its own: the output's metadata is not the input's. + assert ta._as_dataset(amr).layers == 0 + + +def test_average_down_puts_the_fine_values_in_the_coarse_cells_and_leaves_its_input(): + amr = _pulse() + before = [np.asarray(block.cell_data["Gaussian-Pulse"]).copy() + for _, _, block in ta._blocks(amr)] + out = ta.AverageDown(field="Gaussian-Pulse")(amr) + assert isinstance(out, vtk.vtkOverlappingAMR) + for kept, (_, _, block) in zip(before, ta._blocks(amr)): + np.testing.assert_array_equal(np.asarray(block.cell_data["Gaussian-Pulse"]), kept) + expected = _hierarchy(amr) + td.average_down(expected, PULSE) + changed = False + for (_, box), (_, _, block), kept in zip(expected, ta._blocks(out), before): + values = np.asarray(block.cell_data["Gaussian-Pulse"]) + np.testing.assert_array_equal(values, box.fields[PULSE].values.numpy()) + changed |= not np.array_equal(values, kept) + assert changed # the coarse box's refined cells + + +def test_a_field_filter_runs_box_by_box_and_keeps_the_hierarchy(): + """Each block gets the fields the filter computes on it alone; the output is the + input's hierarchy, its blocks shallow copies, and the next filter takes it.""" + amr = _pulse() + out = ta.Elevation()(amr) + assert isinstance(out, vtk.vtkOverlappingAMR) + assert out.GetNumberOfLevels() == amr.GetNumberOfLevels() + for (level, index, block), (_, _, given) in zip(ta._blocks(out), ta._blocks(amr)): + assert out.GetAMRBox(level, index).GetNumberOfCells() == given.number_of_cells + alone = ta.Elevation()(given) + np.testing.assert_array_equal(np.asarray(block.point_data["elevation"]), + np.asarray(alone.point_data["elevation"])) + assert "elevation" not in given.point_data.keys() + # The input's Tack hierarchy with the field: the same topologies, nothing built again. + box = next(box for _, box in ta._as_dataset(out)) + assert "elevation" in box.fields + assert box.topology is next(box for _, box in ta._as_dataset(amr)).topology + table = ta.Statistics(fields=["elevation"])(out) + assert table.GetColumnByName("n").GetValue(0) > 0 + + +def test_a_filter_that_makes_a_dataset_gives_a_partitioned_dataset_per_level(): + """Box by box, the cells finer boxes cover left out: each level's cells are those + the Tack filter keeps of its boxes in the hierarchy, marks and all: the coarse + level's, fewer than its block has. A contour after the points' values is one per + level.""" + amr = _pulse() + hierarchy = _hierarchy(amr) + threshold = ta.Threshold(field=("Gaussian-Pulse", "cells"), lower=-1.0, upper=1.0) + out = threshold(amr) + assert isinstance(out, vtk.vtkPartitionedDataSetCollection) + assert out.GetNumberOfPartitionedDataSets() == amr.GetNumberOfLevels() + for level, pieces in enumerate(hierarchy.levels): + assert out.GetMetaData(level).Get(vtk.vtkCompositeDataSet.NAME()) == f"level{level}" + got = sum(block.number_of_cells for block in out.GetPartitionedDataSet(level)) + expected = sum(td.threshold(box, PULSE, -1.0, 1.0).topology.num_cells + for box in pieces) + assert got == expected + coarse = sum(block.number_of_cells for block in out.GetPartitionedDataSet(0)) + assert 0 < coarse < amr.GetDataSetAsCartesianGrid(0, 0).number_of_cells # refined out + peak = max(np.asarray(block.cell_data["Gaussian-Pulse"]).max() + for _, _, block in ta._blocks(amr)) + surfaces = (ta.Project(field=("Gaussian-Pulse", "cells"), to="points", name="p") + >> ta.Contour(field="p", value=peak / 4))(amr) + assert isinstance(surfaces, vtk.vtkPartitionedDataSetCollection) + assert sum(block.number_of_cells for level in range(amr.GetNumberOfLevels()) + for block in surfaces.GetPartitionedDataSet(level)) > 0 + + +def test_a_hierarchy_is_probed_at_its_points_box_by_box(): + amr = _pulse() + probe = ta.Probe() + probe.SetInputDataObject(0, amr) + probe.SetSourceData(_wavelet((-4, 4, -4, 4, -4, 4))()) + probe.Update() + assert probe.error is None, probe.error + out = probe.GetOutputDataObject(0) + assert isinstance(out, vtk.vtkPartitionedDataSetCollection) + assert all("RTData" in block.point_data.keys() + for level in range(out.GetNumberOfPartitionedDataSets()) + for block in out.GetPartitionedDataSet(level)) + + +def test_filters_whose_answer_box_by_box_would_mislead_refuse_a_hierarchy(): + amr = _pulse() + ids = ta.GenerateGlobalIds() + ids.SetInputDataObject(0, amr) + assert ids.Update() == 0 and isinstance(ids.error, TypeError) + assert "would mislead" in str(ids.error) + lines = ta.Streamlines(field="v", seeds=np.zeros((1, 3))) + lines.SetInputDataObject(0, amr) + assert lines.Update() == 0 and isinstance(lines.error, TypeError) + seeded = ta.Streamlines(field="v", step_size=1e-4, steps=2) + seeded.SetInputDataObject(0, _small()) + seeded.SetSourceData(amr) + assert seeded.Update() == 0 and "port 1" in str(seeded.error) + + +def test_an_amr_hierarchy_is_one_tack_hierarchy_until_a_block_changes(): + """In its shallow copies too; a block's array modified, one added, or a block + replaced, it is converted again (its own modification time follows none).""" + amr = _pulse() + first = ta._as_dataset(amr) + assert ta._as_dataset(amr) is first + copy = amr.NewInstance() + copy.ShallowCopy(amr) + assert ta._as_dataset(copy) is first + amr.GetDataSetAsCartesianGrid(1, 0).cell_data["Gaussian-Pulse"].Modified() + second = ta._as_dataset(amr) + assert second is not first + amr.GetDataSetAsCartesianGrid(1, 0).cell_data["twice"] = 2 * np.asarray( + amr.GetDataSetAsCartesianGrid(1, 0).cell_data["Gaussian-Pulse"]) + third = ta._as_dataset(amr) + assert third is not second + assert "twice" in [box for _, box in third][1].fields # level 1's first box + block = amr.GetDataSetAsCartesianGrid(0, 0).NewInstance() + block.DeepCopy(amr.GetDataSetAsCartesianGrid(0, 0)) + amr.SetDataSet(0, 0, block) + assert ta._as_dataset(amr) is not third + + +def test_an_amr_filter_runs_on_the_device_it_is_given(): + from tack.runtime import dispatch + from tack.runtime.cpu import CPUBackend + + class SecondCPU(CPUBackend): + """A CPU backend posing as another device's, as ``test_partitioned.py`` makes.""" + + display_name = "second device" + + @property + def memory_domain(self): + return ("cpu", 1) + + amr = _pulse() + other = dispatch.as_device(SecondCPU()) + table = ta.Statistics(fields=["Gaussian-Pulse"], device=other)(amr) + once = ta.Statistics(fields=["Gaussian-Pulse"])(amr) + assert all(box.device == other for _, box in ta._as_dataset(amr, other)) + assert table.GetColumnByName("mean").GetValue(0) == pytest.approx( + once.GetColumnByName("mean").GetValue(0), rel=1e-12) + + +# ── Temporal filters: in situ and driving ─────────────────────────── + +TIMES = (0.0, 0.5, 1.0, 1.5, 2.0) + + +def _flow_step(t): + """A step at ``t``: a grid with a velocity ``v`` turning with time and a scalar + ``T`` growing with it.""" + grid = td.uniform_grid((11, 11, 3)) + grid.time = t + x = grid.positions().astype(np.float64) + v = np.zeros_like(x) + v[:, 0], v[:, 1] = 1.0, 0.2 * t + grid.fields["v"] = td.Field(td.H1(grid), tack.array(v.astype(np.float32), tack.vec3())) + grid.fields["T"] = td.Field(td.H1(grid), tack.array((x[:, 0] * (1.0 + t) ** 2) + .astype(np.float32))) + return grid + + +def _pipeline_keys(): + from vtkmodules.vtkCommonExecutionModel import vtkStreamingDemandDrivenPipeline + + return vtkStreamingDemandDrivenPipeline + + +def _steps_source(times=TIMES, no_prior_access=None): + """A source of ``_flow_step``s at ``times`` that records each time it is asked + for (``asked``); with ``no_prior_access``, it says so as an in situ adaptor does + (``NO_PRIOR_TEMPORAL_ACCESS``, that value).""" + from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase + + s = _pipeline_keys() + + class Steps(VTKPythonAlgorithmBase): + def __init__(self): + VTKPythonAlgorithmBase.__init__(self, nInputPorts=0, nOutputPorts=1, + outputType="vtkImageData") + self.asked = [] + self.no_prior_access = no_prior_access + + def RequestInformation(self, request, inInfo, outInfo): + info = outInfo.GetInformationObject(0) + info.Set(s.TIME_STEPS(), list(times), len(times)) + info.Set(s.TIME_RANGE(), [times[0], times[-1]], 2) + info.Set(s.WHOLE_EXTENT(), [0, 10, 0, 10, 0, 2], 6) + if self.no_prior_access is not None: + info.Set(s.NO_PRIOR_TEMPORAL_ACCESS(), self.no_prior_access) + return 1 + + def RequestData(self, request, inInfo, outInfo): + info = outInfo.GetInformationObject(0) + t = info.Get(s.UPDATE_TIME_STEP()) if info.Has(s.UPDATE_TIME_STEP()) else times[0] + self.asked.append(t) + out = self.GetOutputData(outInfo, 0) + out.ShallowCopy(dataset_to_vtk(_flow_step(t))) + out.GetInformation().Set(vtk.vtkDataObject.DATA_TIME_STEP(), t) + return 1 + + return Steps() + + +def _in_situ(make, times=TIMES): + """``make()``'s output after each step handed to it in turn, as in situ.""" + algorithm, outputs = make(), [] + for t in times: + algorithm.SetInputDataObject(0, dataset_to_vtk(_flow_step(t))) + algorithm.Update() + assert algorithm.error is None, algorithm.error + out = algorithm.GetOutputDataObject(0) + copy = out.NewInstance() + copy.DeepCopy(out) + copy.GetInformation().Set(vtk.vtkDataObject.DATA_TIME_STEP(), + out.GetInformation().Get(vtk.vtkDataObject.DATA_TIME_STEP())) + outputs.append(copy) + return outputs + + +def _driven(make, source): + algorithm = make() + algorithm.SetInputConnection(0, source.GetOutputPort()) + return algorithm + + +def _time_of(data): + info = data.GetInformation() + key = vtk.vtkDataObject.DATA_TIME_STEP() + return info.Get(key) if info.Has(key) else None + + +def _same_arrays(a, b, names): + for name in names: + np.testing.assert_array_equal(np.asarray(a.point_data[name]), + np.asarray(b.point_data[name]), err_msg=name) + + +def test_running_statistics_driving_take_every_step_once_and_give_no_time(): + """Driving gives what the steps give in situ after the last, at no time, its output + reporting no time steps; asked again (an update, another time requested) it runs + nothing again.""" + make = lambda: ta.RunningStatistics(fields=["T"]) + source = _steps_source() + stats = _driven(make, source) + stats.Update() + assert stats.error is None, stats.error + out = stats.GetOutputDataObject(0) + assert source.asked == list(TIMES) + assert _time_of(out) is None + s = _pipeline_keys() + assert not stats.GetOutputInformation(0).Has(s.TIME_STEPS()) + assert out.GetFieldData().GetArray("count").GetValue(0) == len(TIMES) + assert sorted(out.point_data.keys()) == ["T_max", "T_mean", "T_min", "T_variance"] + last = _in_situ(make)[-1] + assert _time_of(last) == TIMES[-1] + _same_arrays(out, last, ["T_mean", "T_variance", "T_min", "T_max"]) + expected = np.mean([_flow_step(t).fields["T"].values.numpy() for t in TIMES], axis=0) + np.testing.assert_allclose(np.asarray(out.point_data["T_mean"]), expected, rtol=1e-6) + stats.Update() + stats.UpdateTimeStep(1.0) + assert source.asked == list(TIMES) + stats.time_range = (0.5, 1.5) # a parameter: starts over + stats.Update() + assert source.asked == [*TIMES, 0.5, 1.0, 1.5] + assert stats.GetOutputDataObject(0).GetFieldData().GetArray("count").GetValue(0) == 3 + + +@pytest.mark.parametrize("make", [ + lambda: ta.Pathlines(field="v", seeds=[[1.0, 5.0, 1.0], [2.0, 2.0, 1.0]]), + lambda: ta.Advection(field="v", seeds=[[1.0, 5.0, 1.0], [2.0, 2.0, 1.0]])]) +def test_particles_driving_hold_every_step_and_in_situ_the_steps_so_far(make): + """Driving: pathlines through every step; advection's particles at every step in + one cloud, ``step`` beside ``time``. They are what the same steps give in situ.""" + source = _steps_source() + driven = _driven(make, source) + driven.Update() + assert driven.error is None, driven.error + out = driven.GetOutputDataObject(0) + assert source.asked == list(TIMES) and _time_of(out) is None + in_situ = _in_situ(make) + if isinstance(driven, ta.Pathlines): + last = in_situ[-1] + np.testing.assert_array_equal(np.asarray(out.points.data), np.asarray(last.points.data)) + assert out.number_of_cells == 2 + return + steps = np.asarray(out.point_data["step"]) + assert out.number_of_points == sum(o.number_of_points for o in in_situ) + for k, step in enumerate(in_situ): + at = steps == k + np.testing.assert_array_equal(np.asarray(out.points.data)[at], + np.asarray(step.points.data)) + np.testing.assert_array_equal(np.asarray(out.point_data["time"])[at], + np.full(step.number_of_points, TIMES[k])) + assert _time_of(step) == TIMES[k] + + +def test_streaklines_drive_to_the_time_requested_from_the_steps_taken(): + """At the time requested, stamped with it, as in situ after those steps; a later + time takes only the steps after; an earlier one starts over; the same one + takes none.""" + make = lambda: ta.Streaklines(field="v", seeds=[[1.0, 5.0, 1.0]], releases=2) + source = _steps_source() + streaks = _driven(make, source) + streaks.UpdateTimeStep(1.0) + assert streaks.error is None, streaks.error + assert source.asked == [0.0, 0.5, 1.0] + out = streaks.GetOutputDataObject(0) + assert _time_of(out) == 1.0 + in_situ = _in_situ(make) + np.testing.assert_array_equal(np.asarray(out.points.data), + np.asarray(in_situ[2].points.data)) + streaks.UpdateTimeStep(2.0) + assert source.asked == [0.0, 0.5, 1.0, 1.5, 2.0] + np.testing.assert_array_equal(np.asarray(streaks.GetOutputDataObject(0).points.data), + np.asarray(in_situ[4].points.data)) + streaks.UpdateTimeStep(0.5) + assert source.asked[5:] == [0.0, 0.5] + streaks.UpdateTimeStep(0.5) + assert len(source.asked) == 7 + + +def test_a_time_derivative_takes_the_step_before_once(): + source = _steps_source() + rate = _driven(lambda: ta.TimeDerivative(fields=["T"]), source) + rate.UpdateTimeStep(1.0) + assert rate.error is None, rate.error + assert source.asked == [0.5, 1.0] + out = rate.GetOutputDataObject(0) + assert _time_of(out) == 1.0 and "T" in out.point_data.keys() + t = lambda time: _flow_step(time).fields["T"].values.numpy() + np.testing.assert_allclose(np.asarray(out.point_data["T_dt"]), (t(1.0) - t(0.5)) / 0.5, + rtol=1e-5, atol=1e-5) + rate.UpdateTimeStep(1.5) + assert source.asked == [0.5, 1.0, 1.5] + rate.UpdateTimeStep(0.0) # no step before: no rate + assert source.asked[3:] == [0.0] + assert "T_dt" not in rate.GetOutputDataObject(0).point_data.keys() + in_situ = _in_situ(lambda: ta.TimeDerivative(fields=["T"])) + assert "T_dt" not in in_situ[0].point_data.keys() + rate.UpdateTimeStep(2.0) + _same_arrays(rate.GetOutputDataObject(0), in_situ[4], ["T_dt"]) + + +def test_in_situ_a_step_is_taken_once_and_going_back_or_a_parameter_starts_over(): + stats = ta.RunningStatistics(fields=["T"]) + count = lambda: stats.GetOutputDataObject(0).GetFieldData().GetArray("count").GetValue(0) + for t in (0.0, 0.5, 0.5): # the same step again: not taken + stats.SetInputDataObject(0, dataset_to_vtk(_flow_step(t))) + stats.Update() + assert stats.error is None and count() == 2 + assert _time_of(stats.GetOutputDataObject(0)) == 0.5 + stats.SetInputDataObject(0, dataset_to_vtk(_flow_step(0.0))) + stats.Update() + assert count() == 1 # back in time: started over + stats.SetInputDataObject(0, dataset_to_vtk(_flow_step(0.5))) + stats.Update() + assert count() == 2 + stats.fields = [("T", "points")] + stats.SetInputDataObject(0, dataset_to_vtk(_flow_step(1.0))) + stats.Update() + assert count() == 1 # a parameter changed + + +def test_an_in_situ_adaptor_keeps_a_filter_from_driving(): + """With VTK's NO_PRIOR_TEMPORAL_ACCESS on its input a filter takes the step it is + given, as does mode='in situ' anywhere; RESET starts it over.""" + s = _pipeline_keys() + source = _steps_source(no_prior_access=s.NO_PRIOR_TEMPORAL_ACCESS_CONTINUE) + stats = _driven(lambda: ta.RunningStatistics(fields=["T"]), source) + for t in (0.0, 0.5, 1.0): + stats.UpdateTimeStep(t) + assert stats.error is None, stats.error + assert source.asked == [0.0, 0.5, 1.0] + out = stats.GetOutputDataObject(0) + assert out.GetFieldData().GetArray("count").GetValue(0) == 3 and _time_of(out) == 1.0 + source.no_prior_access = s.NO_PRIOR_TEMPORAL_ACCESS_RESET + source.Modified() + stats.UpdateTimeStep(1.5) + assert stats.GetOutputDataObject(0).GetFieldData().GetArray("count").GetValue(0) == 1 + plain = _steps_source() + forced = _driven(lambda: ta.RunningStatistics(fields=["T"], mode="in situ"), plain) + forced.UpdateTimeStep(1.0) + assert plain.asked == [1.0] + + +def test_seeds_on_port_1_and_trees_refused(): + seeds = vtk.vtkPolyData() + seeds.points = vtk.vtkPoints(data=np.array([[1.0, 5.0, 1.0]])) + source = _steps_source() + paths = ta.Pathlines(field="v") + paths.SetInputConnection(0, source.GetOutputPort()) + paths.SetSourceData(seeds) + paths.Update() + assert paths.error is None, paths.error + assert paths.GetOutputDataObject(0).number_of_cells == 1 + tree = vtk.vtkMultiBlockDataSet() + tree.SetNumberOfBlocks(1) + tree.SetBlock(0, dataset_to_vtk(_flow_step(0.0))) + stats = ta.RunningStatistics(fields=["T"]) + stats.SetInputDataObject(0, tree) + assert stats.Update() == 0 and isinstance(stats.error, TypeError)