diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 2d624617f..cf99a1c4d 100644 --- a/doc/docs/Python_User_Interface.md +++ b/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). diff --git a/python/simulation.py b/python/simulation.py index cf8b1d3fd..c81c68c58 100644 --- a/python/simulation.py +++ b/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) diff --git a/python/tests/test_simulation.py b/python/tests/test_simulation.py index 5864f7d31..34340ecb7 100644 --- a/python/tests/test_simulation.py +++ b/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) diff --git a/python/tests/test_source.py b/python/tests/test_source.py index acadf8815..c4505d446 100644 --- a/python/tests/test_source.py +++ b/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): diff --git a/src/vec.cpp b/src/vec.cpp index 920bcf912..b5db18574 100644 --- a/src/vec.cpp +++ b/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) {