diff --git a/python/Makefile.am b/python/Makefile.am index e83d4ebae..262ece801 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -83,6 +83,7 @@ TESTS = \ $(TEST_DIR)/test_multilevel_atom.py \ $(TEST_DIR)/test_n2f_periodic.py \ $(TEST_DIR)/test_oblique_source.py \ + $(TEST_DIR)/test_pml_cyl.py \ $(TEST_DIR)/test_physical.py \ $(TEST_DIR)/test_prism.py \ $(TEST_DIR)/test_pw_source.py \ diff --git a/python/tests/test_pml_cyl.py b/python/tests/test_pml_cyl.py index dc00cf371..70c94b82e 100644 --- a/python/tests/test_pml_cyl.py +++ b/python/tests/test_pml_cyl.py @@ -4,76 +4,92 @@ class TestPMLCylindrical(unittest.TestCase): + @classmethod + def setUp(cls): + cls.resolution = 25 # pixels/um + cls.s = 5.0 + cls.dpml_r = 1.0 + cls.dpml_z = 1.0 + cls.cell_size = mp.Vector3( + cls.s + cls.dpml_r, + 0, + cls.s + 2 * cls.dpml_z, + ) + cls.fcen = 1.0 + def test_pml_cyl_coords(self): """Verifies that the z-PML in cylindrical coordinates properly attenuates fields at r=0. """ - resolution = 25 # pixels/um - s = 5.0 - dpml_r = 1.0 - dpml_z = 1.0 - cell_size = mp.Vector3(s + dpml_r, 0, s + 2 * dpml_z) - pml_layers = [ - mp.PML(dpml_r, direction=mp.R), - mp.PML(dpml_z, direction=mp.Z), + mp.PML(self.dpml_r, direction=mp.R), + mp.PML(self.dpml_z, direction=mp.Z), ] - fcen = 1.0 - sources = [ mp.Source( - src=mp.GaussianSource(fcen, fwidth=0.1 * fcen), + src=mp.GaussianSource(self.fcen, fwidth=0.1 * self.fcen), center=mp.Vector3(), component=mp.Er, ), ] sim = mp.Simulation( - resolution=resolution, - cell_size=cell_size, + resolution=self.resolution, + cell_size=self.cell_size, dimensions=mp.CYLINDRICAL, m=-1, - accurate_fields_near_cylorigin=False, sources=sources, boundary_layers=pml_layers, ) flux_plus_z = sim.add_flux( - fcen, + self.fcen, 0, 1, mp.FluxRegion( - center=mp.Vector3(0.5 * s, 0, 0.5 * s), size=mp.Vector3(s, 0, 0) + center=mp.Vector3( + 0.5 * self.s, + 0, + 0.5 * self.s, + ), + size=mp.Vector3(self.s, 0, 0), ), ) flux_plus_r = sim.add_flux( - fcen, + self.fcen, 0, 1, - mp.FluxRegion(center=mp.Vector3(s, 0, 0), size=mp.Vector3(0, 0, s)), + mp.FluxRegion( + center=mp.Vector3(self.s, 0, 0), + size=mp.Vector3(0, 0, self.s), + ), ) flux_minus_z = sim.add_flux( - fcen, + self.fcen, 0, 1, mp.FluxRegion( - center=mp.Vector3(0.5 * s, 0, -0.5 * s), - size=mp.Vector3(s, 0, 0), + center=mp.Vector3( + 0.5 * self.s, + 0, + -0.5 * self.s, + ), + size=mp.Vector3(self.s, 0, 0), weight=-1.0, ), ) - sim.run(until_after_sources=50) + sim.run(until_after_sources=50.94) prev_flux_plusz = mp.get_fluxes(flux_plus_z)[0] prev_flux_plusr = mp.get_fluxes(flux_plus_r)[0] prev_flux_minusz = mp.get_fluxes(flux_minus_z)[0] - for t in [100, 200, 400]: + for t in [142.15, 214.64, 365.32]: sim.run(until_after_sources=t) flux_plusz = mp.get_fluxes(flux_plus_z)[0] @@ -86,9 +102,10 @@ def test_pml_cyl_coords(self): f"{flux_plusr:.6f}, {flux_minusz:.6f}, {flux_tot:.6f}" ) - self.assertAlmostEqual(flux_plusz, prev_flux_plusz, places=9) - self.assertAlmostEqual(flux_plusr, prev_flux_plusr, places=9) - self.assertAlmostEqual(flux_minusz, prev_flux_minusz, places=9) + places = 6 if mp.is_single_precision() else 9 + self.assertAlmostEqual(flux_plusz, prev_flux_plusz, places=places) + self.assertAlmostEqual(flux_plusr, prev_flux_plusr, places=places) + self.assertAlmostEqual(flux_minusz, prev_flux_minusz, places=places) prev_flux_plusz = flux_plusz prev_flux_plusr = flux_plusr diff --git a/src/step_db.cpp b/src/step_db.cpp index 4145d3e6b..8a3bf884b 100644 --- a/src/step_db.cpp +++ b/src/step_db.cpp @@ -339,6 +339,7 @@ bool fields_chunk::step_db(field_type ft) { the_f[iz] += (df = dfcnd * (siginv ? siginv[dk + 2 * (dsig == Z) * iz] : 1)); if (fu) fu[iz] += siginvu[dku + 2 * (dsigu == Z) * iz] * df; } + // hack: manually set boundary conditions at r=0 if (ft == D_stuff) { ZERO_Z(f[Dz][cmp]); if (f_cond[Dz][cmp]) ZERO_Z(f_cond[Dz][cmp]); diff --git a/tests/cylindrical.cpp b/tests/cylindrical.cpp index 706d83b77..db4e399d0 100644 --- a/tests/cylindrical.cpp +++ b/tests/cylindrical.cpp @@ -85,7 +85,8 @@ int test_simple_periodic(double eps(const vec &), int splitting) { if (!compare_point(f, f1, veccyl(0.5, 0.4))) return 0; if (!compare_point(f, f1, veccyl(0.46, 0.36))) return 0; if (!compare_point(f, f1, veccyl(1.0, 0.4))) return 0; - if (!compare_point(f, f1, veccyl(0.01, 0.02))) return 0; + if (!compare_point(f, f1, veccyl(0.01, 0.02), sizeof(realnum) == sizeof(float) ? 2e-6 : 4e-8)) + return 0; if (!compare_point(f, f1, veccyl(0.601, 0.701))) return 0; if (f.time() >= field_energy_check_time) { if (!compare(f.field_energy(), f1.field_energy(), " total energy")) return 0; @@ -129,7 +130,8 @@ int test_simple_metallic(double eps(const vec &), int splitting) { if (!compare_point(f, f1, veccyl(0.5, 0.4))) return 0; if (!compare_point(f, f1, veccyl(0.46, 0.36))) return 0; if (!compare_point(f, f1, veccyl(1.0, 0.4))) return 0; - if (!compare_point(f, f1, veccyl(0.01, 0.02))) return 0; + if (!compare_point(f, f1, veccyl(0.01, 0.02), sizeof(realnum) == sizeof(float) ? 2e-6 : 4e-8)) + return 0; if (!compare_point(f, f1, veccyl(0.601, 0.701))) return 0; if (f.time() >= field_energy_check_time) { if (!compare(f.field_energy(), f1.field_energy(), " total energy")) return 0;