increase tolerance of check for r≈0 fields of cylindrical.cpp for sin…

…gle-precision builds
NanoComp · Jan 26, 2023 · f4c26a4 · f4c26a4
1 parent 523bb9d
commit f4c26a4
Show file tree

Hide file tree

Showing 4 changed files with 49 additions and 28 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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;