From 50cee9154ce372fcce66cd928e9959d67230f103 Mon Sep 17 00:00:00 2001 From: Max Date: Wed, 28 Oct 2026 16:21:00 +0100 Subject: [PATCH 1/2] Add the current deposit along paths to the edge space. --- CHANGELOG.md | 5 ++ src/patchwork/edges.py | 86 ++++++++++++++++++++++++++ src/patchwork/tests/unit/test_edges.py | 71 +++++++++++++++++++++ tutorials/07_edge_space.ipynb | 35 ++++++++++- 4 files changed, 195 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f0b58b7..810b468 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,11 @@ format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). for `G phi`) and `deposit(points, vectors)`, its transpose. `CompositeNodes.point_cells` gives the leaf cell of points and where in it they are. Tutorial: the edge space. +- `CompositeEdges.deposit_paths(start, end, charges)`: the current of charges moving along + straight paths, each path split at every cell face it crosses (patch edges and periodic + seams included) and integrated exactly against the edge functions. The charge is conserved + to round-off: `(G T)^T j` is the change of the deposited nodal charge. A NumPy reference. + - `AMRHierarchy.wrap(points)` and `wrapped(points)`: points beyond a periodic face moved back into the domain, in place or as a copy (the array itself when nothing needs wrapping), and `AMRHierarchy.node_grid(patch)`: a patch's nodes as an `mpifield.BlockGrid`, for compiled diff --git a/src/patchwork/edges.py b/src/patchwork/edges.py index 7ef5a1b..5e2441f 100644 --- a/src/patchwork/edges.py +++ b/src/patchwork/edges.py @@ -273,3 +273,89 @@ def deposit(self, points: np.ndarray, vectors: np.ndarray) -> np.ndarray: vectors = np.asarray(vectors, dtype=float) x, y = self.evaluation(points) return x.T @ vectors[0] + y.T @ vectors[1] + + def _pieces( + self, start: np.ndarray, end: np.ndarray + ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """The straight paths from ``start`` to ``end`` ``(2, k)``, cut at the finest lattice. + + Returns the path of every piece and its two ends ``(2, pieces)``. A + piece lies in one cell of the finest lattice, and so in one leaf cell. + """ + nodes = self.nodes + origin = nodes._origin[:, None] + spacing = nodes._spacing[:, None] + lower = (start - origin) / spacing + upper = (end - origin) / spacing + count = start.shape[1] + paths, times = [np.arange(count)] * 2, [np.zeros(count), np.ones(count)] + for axis in range(2): + a, b = lower[axis], upper[axis] + first = np.floor(np.minimum(a, b)).astype(np.int64) + 1 + last = np.ceil(np.maximum(a, b)).astype(np.int64) - 1 + crossings = np.maximum(last - first + 1, 0) + path = np.repeat(np.arange(count), crossings) + offset = np.repeat(np.cumsum(crossings) - crossings, crossings) + line = first[path] + np.arange(path.size) - offset + paths.append(path) + times.append((line - a[path]) / (b - a)[path]) + path = np.concatenate(paths) + time = np.concatenate(times) + order = np.lexsort((time, path)) + path, time = path[order], time[order] + # consecutive times of one path bound a piece + same = path[1:] == path[:-1] + path, begin, finish = path[1:][same], time[:-1][same], time[1:][same] + step = end - start + return ( + path, + start[:, path] + begin * step[:, path], + start[:, path] + finish * step[:, path], + ) + + def deposit_paths( + self, start: np.ndarray, end: np.ndarray, charges: np.ndarray + ) -> np.ndarray: + """The current of charges moving along straight paths, integrated over the paths: ``(n1,)``. + + Every charge goes from ``start`` to ``end`` (``(2, k)``, not wrapped into + the domain on a periodic axis); its path is split at every cell face it + crosses, patch edges included, and each piece is integrated exactly + against the edge functions of its leaf cell: ``q * integral of w . dl``. + Then the charge is conserved exactly: ``(G T)^T j`` is the change of the + charge on the nodes, ``T^T`` of :meth:`CompositeNodes.deposit + ` at the ends minus at the + starts, and ``j @ e`` is the work of the field ``e`` on the charges. The + current density times the time step, summed over a step's paths, is + ``j``. Pieces outside the domain deposit nothing. A NumPy reference, + one pass over all pieces. + """ + start = np.asarray(start, dtype=float) + end = np.asarray(end, dtype=float) + charges = np.asarray(charges, dtype=float).reshape(-1) + if start.shape != end.shape or start.shape != (2, charges.size): + raise ValueError( + f"paths from {start.shape} to {end.shape} for {charges.size} charges" + ) + path, a, b = self._pieces(start, end) + nodes = self.nodes + # the edge functions are linear in a leaf cell: the midpoint integrates them + leaf, fraction, inside = nodes.point_cells(0.5 * (a + b)) + size = nodes.cell_size[leaf] + charge = np.where(inside, charges[path], 0.0) + step = b - a + current = np.zeros(self.num_edges) + for direction, (low, high) in ((0, (0, 1)), (1, (2, 3))): + across = fraction[1 - direction] + along = charge * step[direction] / size[:, direction] + current += np.bincount( + self.cell_edges[leaf, low], + weights=(1.0 - across) * along, + minlength=self.num_edges, + ) + current += np.bincount( + self.cell_edges[leaf, high], + weights=across * along, + minlength=self.num_edges, + ) + return current diff --git a/src/patchwork/tests/unit/test_edges.py b/src/patchwork/tests/unit/test_edges.py index 9dd297e..eacb11f 100644 --- a/src/patchwork/tests/unit/test_edges.py +++ b/src/patchwork/tests/unit/test_edges.py @@ -163,3 +163,74 @@ def test_the_edge_space_is_two_dimensional() -> None: nodes = CompositeNodes(AMRHierarchy(base_shape=(8,))) with pytest.raises(NotImplementedError, match="two-dimensional"): CompositeEdges(nodes) + + +def _continuity_error(edges: CompositeEdges, start, end, charges) -> float: + """``(G T)^T j`` against the change of the composite nodal charge, relative.""" + nodes = edges.nodes + current = edges.deposit_paths(start, end, charges) + divergence = (edges.gradient @ nodes.reduction).T @ current + change = nodes.reduction.T @ ( + nodes.deposit(end, charges) - nodes.deposit(start, charges) + ) + return float(np.abs(divergence - change).max() / np.abs(charges).max()) + + +@pytest.mark.parametrize( + ("start", "end"), + [ + ((0.03, 0.03), (0.05, 0.06)), # inside a coarse cell + ((0.03, 0.03), (0.2, 0.05)), # across faces of coarse cells + ((0.45, 0.4), (0.55, 0.62)), # across the edge of the level-1 patch + ((0.4, 0.4), (1.3, 0.85)), # through both patches and out + ((1.0, 0.97), (1.1, 1.08)), # across the periodic seam + ((0.6, 0.5), (0.6, 0.7)), # along a lattice line + ((0.5, 0.5), (0.75, 0.75)), # through lattice corners + ], +) +def test_the_path_deposit_conserves_the_charge(start, end) -> None: + edges = _edges("nested", periodic=(False, True)) + start, end = np.array(start)[:, None], np.array(end)[:, None] + assert _continuity_error(edges, start, end, np.array([1.5])) < 1e-12 + + +@pytest.mark.parametrize("layout", LAYOUTS) +def test_many_random_paths_conserve_the_charge(layout: str) -> None: + edges = _edges(layout, periodic=(False, True)) + rng = np.random.default_rng(4) + start = _points(500, seed=9) * np.array([[0.9], [1.0]]) + np.array([[0.1], [0.0]]) + end = start + rng.normal(scale=0.15, size=start.shape) + end[0] = np.clip(end[0], 0.0, 2.0) # bounded in x, periodic in y + charges = rng.normal(size=500) + assert _continuity_error(edges, start, end, charges) < 1e-12 + + +def test_the_work_of_a_gradient_field_is_the_potential_difference() -> None: + edges = _edges("level jump", periodic=(False, True)) + nodes = edges.nodes + phi = _potential(nodes) + start = np.array([[0.3, 0.9, 1.6], [0.2, 0.95, 0.4]]) + end = np.array([[1.4, 1.1, 0.7], [0.65, 1.3, 0.55]]) + charges = np.array([1.0, -2.0, 0.5]) + work = edges.deposit_paths(start, end, charges) @ (edges.gradient @ phi) + expected = charges @ (nodes.evaluation(end) @ phi - nodes.evaluation(start) @ phi) + assert work == pytest.approx(expected, rel=1e-12) + + +def test_a_uniform_field_does_the_work_of_its_displacement() -> None: + edges = _edges("two patches") + field = np.array([2.0, -0.5]) + e = edges.length * field[edges.direction] + start = np.array([[0.1, 1.9], [0.1, 0.8]]) + end = np.array([[1.7, 0.2], [0.9, 0.15]]) + charges = np.array([0.7, 1.1]) + work = edges.deposit_paths(start, end, charges) @ e + assert work == pytest.approx(float(charges @ (field @ (end - start))), rel=1e-12) + + +def test_a_path_that_stays_put_deposits_nothing_and_shapes_are_checked() -> None: + edges = _edges("nested") + point = np.array([[0.4], [0.3]]) + np.testing.assert_array_equal(edges.deposit_paths(point, point, np.ones(1)), 0.0) + with pytest.raises(ValueError, match="paths"): + edges.deposit_paths(point, np.zeros((2, 2)), np.ones(1)) diff --git a/tutorials/07_edge_space.ipynb b/tutorials/07_edge_space.ipynb index 7b19af4..def4da5 100644 --- a/tutorials/07_edge_space.ipynb +++ b/tutorials/07_edge_space.ipynb @@ -145,8 +145,7 @@ "\n", "`edges.deposit(points, vectors)` puts vectors at points on the edges with the weights of the\n", "gather, whose transpose it is: the work of the gathered field on the markers' currents is the\n", - "current's work on the edges. Charge is conserved exactly when the current is integrated along\n", - "each marker's path; that is the next step of the scheme." + "current's work on the edges." ] }, { @@ -162,6 +161,38 @@ "gathered = np.array([matrix @ e for matrix in edges.evaluation(points)])\n", "print(edges.deposit(points, current) @ e, np.sum(current * gathered))" ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## The current along paths\n", + "\n", + "`edges.deposit_paths(start, end, charges)` integrates the current of charges moving along\n", + "straight paths against the edge functions, each path split at every cell face it crosses,\n", + "patch edges included. The charge is then conserved exactly: the divergence of the current,\n", + "`(G T)^T j`, is the change of the charge deposited on the nodes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "start = np.vstack((rng.random(500), rng.random(500)))\n", + "end = start + rng.normal(scale=0.1, size=start.shape)\n", + "end[0] = np.clip(end[0], 0.0, 1.0) # bounded in x, periodic in y\n", + "charges = rng.normal(size=500)\n", + "current = edges.deposit_paths(start, end, charges)\n", + "divergence = gradient.T @ current\n", + "change = nodes.reduction.T @ (\n", + " nodes.deposit(end, charges) - nodes.deposit(start, charges)\n", + ")\n", + "print(\"largest difference:\", np.abs(divergence - change).max())" + ] } ], "metadata": { From a47a685caab86aa7dff454061f5b7a48f32ca097 Mon Sep 17 00:00:00 2001 From: Max Date: Wed, 28 Oct 2026 17:04:25 +0100 Subject: [PATCH 2/2] Run the MPI tests in CI without UCX, which fails on some runners. --- .github/workflows/test_pytest.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/test_pytest.yml b/.github/workflows/test_pytest.yml index d3a0f67..a58db26 100644 --- a/.github/workflows/test_pytest.yml +++ b/.github/workflows/test_pytest.yml @@ -76,6 +76,11 @@ jobs: - name: Run the PETSc and MPI tests env: OMP_NUM_THREADS: "1" + # Open MPI's UCX cannot open the network card of some runners (Azure MANA: + # "uct_iface_open ... Address not valid") and the ranks then fail or hang; + # on one machine the shared-memory transport of ob1 is all the tests need + OMPI_MCA_pml: ob1 + OMPI_MCA_osc: ^ucx run: | python -c "from petsc4py import PETSc; print(PETSc.Sys.getVersion())" python -m pytest -ra src/patchwork/tests/unit/test_distributed.py src/patchwork/tests/unit/test_mpi.py