Add the current deposit along paths to the edge space

.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

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

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
+        <patchwork.nodal.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

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))

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": {