[2/13][Adjoint Module] Round the grid origin to a whole pixel

doc/docs/Python_User_Interface.md

@@ -141,7 +141,9 @@ Python. `Vector3` is a `meep` class.
   the coordinate system used in Meep (for example, to put the origin at the
   corner).  Passing `geometry_center=c` is equivalent to adding the `c` vector to
   the coordinates of every other object in the simulation, i.e. `c` becomes the
-  new origin that other objects are defined with respect to.
+  new origin that other objects are defined with respect to. `c` is rounded to the
+  nearest pixel (a multiple of `1/resolution`) so that the grid is not shifted
+  relative to the Yee lattice.
 
 + **`sources` [ list of `Source` class ]** — Specifies the current sources to be
   present in the simulation. Defaults to none (empty list).

python/simulation.py

@@ -1350,7 +1350,9 @@ def __init__(
           the coordinate system used in Meep (for example, to put the origin at the
           corner).  Passing `geometry_center=c` is equivalent to adding the `c` vector to
           the coordinates of every other object in the simulation, i.e. `c` becomes the
-          new origin that other objects are defined with respect to.
+          new origin that other objects are defined with respect to. `c` is rounded to the
+          nearest pixel (a multiple of `1/resolution`) so that the grid is not shifted
+          relative to the Yee lattice.
 
         + **`sources` [ list of `Source` class ]** — Specifies the current sources to be
           present in the simulation. Defaults to none (empty list).
@@ -1788,6 +1790,14 @@ def _create_grid_volume(self, k: Vector3Type = None):
         else:
             raise ValueError(f"Unsupported dimentionality: {dims}")
 
+        # the grid origin is rounded to a whole pixel (see grid_volume::set_origin),
+        # so round geometry_center the same way to keep it consistent with the grid
+        self.geometry_center = Vector3(
+            *(
+                round(c * self.resolution) / self.resolution
+                for c in self.geometry_center
+            )
+        )
         gv.center_origin()
         gv.shift_origin(
             py_v3_to_vec(self.dimensions, self.geometry_center, self.is_cylindrical)

python/tests/test_simulation.py

@@ -715,6 +715,46 @@ def print_field(sim):
 
         self.assertAlmostEqual(result[0], -0.0599602798684155)
 
+    def test_geometry_center_rounded_to_pixel(self):
+        # a geometry_center that is not a whole number of pixels used to give the grid an
+        # odd origin (in half-pixel units), which put sources and monitors half a pixel
+        # away from the fields and materials and made results depend on the chunking
+        resolution = 8
+
+        def get_ez(geometry_center):
+            sim = mp.Simulation(
+                resolution=resolution,
+                cell_size=mp.Vector3(4, 4),
+                geometry_center=geometry_center,
+                geometry=[
+                    mp.Block(
+                        center=mp.Vector3(0.59, 0.21),
+                        size=mp.Vector3(0.8, 1.3),
+                        material=mp.Medium(epsilon=6),
+                    )
+                ],
+                sources=[
+                    mp.Source(
+                        mp.GaussianSource(1.0, fwidth=0.5),
+                        component=mp.Ez,
+                        center=mp.Vector3(0.0234375, 0.0078125),
+                    )
+                ],
+                eps_averaging=False,
+            )
+            sim.run(until=5)
+            io = sim.fields.gv.little_corner()
+            self.assertEqual((io.x() % 2, io.y() % 2), (0, 0))
+            return (
+                sim.get_field_point(mp.Ez, mp.Vector3(1.3, -0.4)),
+                sim.geometry_center,
+            )
+
+        ez, geometry_center = get_ez(mp.Vector3(0.0625, 0.3))
+        ez_rounded, _ = get_ez(mp.Vector3(0, 0.25))
+        self.assertTrue(geometry_center.close(mp.Vector3(0, 0.25)))
+        self.assertEqual(ez, ez_rounded)
+
     def test_timing_data(self):
         resolution = 20
         cell_size = mp.Vector3(10, 10)

python/tests/test_source.py

@@ -125,6 +125,114 @@ def test_gaussian_fourier_transform_continuous(self):
             self.assertEqual(src.fourier_transform(f), src.swigobj.fourier_transform(f))
 
 
+class TestSourceChunkBoundary(unittest.TestCase):
+    @staticmethod
+    def _get_source(
+        cell_size,
+        geometry_center,
+        source_center,
+        source_size,
+        chunk_layout,
+        run_until=None,
+    ):
+        sim = mp.Simulation(
+            cell_size=cell_size,
+            geometry_center=geometry_center,
+            resolution=8,
+            chunk_layout=chunk_layout,
+            sources=[
+                mp.Source(
+                    mp.ContinuousSource(1.0),
+                    component=mp.Ez,
+                    center=source_center,
+                    size=source_size,
+                )
+            ],
+        )
+        try:
+            sim.init_sim()
+            source = sim.get_source(mp.Ez)
+            energy = None
+            if run_until is not None:
+                sim.run(until=run_until)
+                energy = sim.field_energy_in_box(
+                    mp.Volume(center=geometry_center, size=cell_size)
+                )
+            return source, energy
+        finally:
+            sim.reset_meep()
+
+    def test_zero_thickness_source_at_chunk_boundary(self):
+        cell_size = mp.Vector3(2.5, 2.25, 2.75)
+        geometry_center = mp.Vector3(0.0625, -0.0625, 0.0625)
+        one_chunk = mp.BinaryPartition(data=0)
+        split_at_source = mp.BinaryPartition(data=[(mp.Z, 0.0), 0, 1])
+
+        # The off-boundary case is a control that the split alone does not
+        # change deposition. The boundary case exercises interpolation onto
+        # both Ez planes adjacent to the split.
+        for source_z in (0.25, 0.0):
+            with self.subTest(source_z=source_z):
+                source_center = mp.Vector3(0.0625, -0.0625, source_z)
+                source_size = mp.Vector3(1.5, 1.25, 0)
+                expected_source, expected_energy = self._get_source(
+                    cell_size,
+                    geometry_center,
+                    source_center,
+                    source_size,
+                    one_chunk,
+                    run_until=0.25 if source_z == 0 else None,
+                )
+                actual_source, actual_energy = self._get_source(
+                    cell_size,
+                    geometry_center,
+                    source_center,
+                    source_size,
+                    split_at_source,
+                    run_until=0.25 if source_z == 0 else None,
+                )
+
+                self.assertTrue(np.any(expected_source != 0))
+                np.testing.assert_array_equal(actual_source != 0, expected_source != 0)
+                np.testing.assert_array_equal(actual_source, expected_source)
+                if expected_energy is not None:
+                    # Source introspection can include non-owned storage that
+                    # never injects into the evolved fields.
+                    rtol = 1e-5 if mp.is_single_precision() else 1e-12
+                    np.testing.assert_allclose(
+                        actual_energy, expected_energy, rtol=rtol, atol=0
+                    )
+
+    def test_point_source_at_chunk_boundary_2d(self):
+        args = (
+            mp.Vector3(2.5, 2.25),
+            mp.Vector3(0.0625, -0.0625),
+            mp.Vector3(0.0625, -0.0625),
+            mp.Vector3(),
+        )
+        expected, _ = self._get_source(*args, mp.BinaryPartition(data=0))
+        actual, _ = self._get_source(
+            *args, mp.BinaryPartition(data=[(mp.X, 0.0), 0, 1])
+        )
+
+        self.assertTrue(np.any(expected != 0))
+        np.testing.assert_array_equal(actual != 0, expected != 0)
+        np.testing.assert_array_equal(actual, expected)
+
+    def test_point_source_on_owned_chunk_boundary_2d(self):
+        # This centered-grid control guards the historical failure mode where
+        # both chunks deposited the same already-owned source sample.
+        args = (mp.Vector3(2, 2), mp.Vector3(), mp.Vector3(), mp.Vector3())
+        expected, _ = self._get_source(*args, mp.BinaryPartition(data=0))
+        actual, _ = self._get_source(
+            *args, mp.BinaryPartition(data=[(mp.X, 0.0), 0, 1])
+        )
+
+        self.assertTrue(np.any(expected != 0))
+        np.testing.assert_array_equal(actual != 0, expected != 0)
+        np.testing.assert_array_equal(actual, expected)
+
+
 class TestSourceTypemaps(unittest.TestCase):
     def setUp(self):
         def dummy_eps(v):

src/vec.cpp

@@ -45,7 +45,12 @@ void grid_volume::set_origin(direction d, int o) {
   origin = operator[](io); // adjust origin to match io
 }
 
-void grid_volume::set_origin(const vec &o) { set_origin(round_vec(o)); }
+void grid_volume::set_origin(const vec &o) {
+  // round to a whole pixel: the Yee lattice assumes io is even (see icenter)
+  ivec io2(dim);
+  LOOP_OVER_DIRECTIONS(dim, d) { io2.set_direction(d, 2 * my_round(o.in_direction(d) * a)); }
+  set_origin(io2);
+}
 
 const char *dimension_name(ndim dim) {
   switch (dim) {