Fix bug for z-PML in cylindrical coordinates for m = 0, ±1 (#2383)

* zero the Br fields at r=0 for m=±1

* increase tolerance of check for r≈0 fields of cylindrical.cpp for single-precision builds

* revert changes in step_db.cpp and remove special case of update E from D from r=0 in update_eh.cpp

* slightly increase tolerance of gradient check in test_adjoint_cyl.py

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              \

python/tests/test_adjoint_cyl.py

@@ -201,7 +201,7 @@ def test_adjoint_solver_cyl_n2f_fields(self, m, far_x):
         adj_scale = (dp[None, :] @ adjsol_grad).flatten()
         fd_grad = S12_perturbed - S12_unperturbed
         print(f"Directional derivative -- adjoint solver: {adj_scale}, FD: {fd_grad}")
-        tol = 0.2 if mp.is_single_precision() else 0.1
+        tol = 0.3
         self.assertClose(adj_scale, fd_grad, epsilon=tol)
 
 

python/tests/test_pml_cyl.py

@@ -0,0 +1,118 @@
+import unittest
+import meep as mp
+import numpy as np
+
+
+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.
+        """
+
+        pml_layers = [
+            mp.PML(self.dpml_r, direction=mp.R),
+            mp.PML(self.dpml_z, direction=mp.Z),
+        ]
+
+        sources = [
+            mp.Source(
+                src=mp.GaussianSource(self.fcen, fwidth=0.1 * self.fcen),
+                center=mp.Vector3(),
+                component=mp.Er,
+            ),
+        ]
+
+        sim = mp.Simulation(
+            resolution=self.resolution,
+            cell_size=self.cell_size,
+            dimensions=mp.CYLINDRICAL,
+            m=-1,
+            sources=sources,
+            boundary_layers=pml_layers,
+        )
+
+        flux_plus_z = sim.add_flux(
+            self.fcen,
+            0,
+            1,
+            mp.FluxRegion(
+                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(
+            self.fcen,
+            0,
+            1,
+            mp.FluxRegion(
+                center=mp.Vector3(self.s, 0, 0),
+                size=mp.Vector3(0, 0, self.s),
+            ),
+        )
+
+        flux_minus_z = sim.add_flux(
+            self.fcen,
+            0,
+            1,
+            mp.FluxRegion(
+                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.94)
+
+        prev_flux = [
+            mp.get_fluxes(flux_plus_z)[0],
+            mp.get_fluxes(flux_plus_r)[0],
+            mp.get_fluxes(flux_minus_z)[0],
+        ]
+
+        for t in [142.15, 214.64, 365.32]:
+            sim.run(until_after_sources=t)
+
+            cur_flux = [
+                mp.get_fluxes(flux_plus_z)[0],
+                mp.get_fluxes(flux_plus_r)[0],
+                mp.get_fluxes(flux_minus_z)[0],
+            ]
+            cur_flux_str = ", ".join(f"{c:.6f}" for c in cur_flux)
+            flux_tot = sum(cur_flux)
+
+            print(f"flux:, {sim.meep_time()}, {cur_flux_str}, {flux_tot:.6f}")
+
+            places = 6 if mp.is_single_precision() else 9
+            for i in range(len(cur_flux)):
+                self.assertAlmostEqual(
+                    prev_flux[i],
+                    cur_flux[i],
+                    places=places,
+                )
+                prev_flux[i] = cur_flux[i]
+
+
+if __name__ == "__main__":
+    unittest.main()

src/update_eh.cpp

@@ -197,34 +197,6 @@ bool fields_chunk::update_eh(field_type ft, bool skip_w_components) {
     }
   }
 
-  /* Do annoying special cases for r=0 in cylindrical coords.  Note
-     that this only really matters for field output; the Ez and Ep
-     components at r=0 don't usually affect the fields elsewhere
-     because of the form of Maxwell's equations in cylindrical coords. */
-  // (FIXME: handle Kerr case?  Do we care about auxiliary PML fields here?)
-  if (gv.dim == Dcyl && gv.origin_r() == 0.0) DOCMP FOR_FT_COMPONENTS(ft, ec) {
-      if (f[ec][cmp] && (ec == Ep || ec == Ez || ec == Hr)) {
-        component dc = field_type_component(ft2, ec);
-        if (f[ec][cmp] == f[dc][cmp]) continue;
-        const int yee_idx = gv.yee_index(ec);
-        const int d_ec = component_direction(ec);
-        const int sR = gv.stride(R), nZ = gv.num_direction(Z);
-        realnum *E = f[ec][cmp];
-        const realnum *D = f_minus_p[dc][cmp] ? f_minus_p[dc][cmp] : f[dc][cmp];
-        const realnum *chi1inv = s->chi1inv[ec][d_ec];
-        if (chi1inv)
-          for (int iZ = 0; iZ < nZ; iZ++) {
-            const int i = yee_idx + iZ - sR;
-            E[i] = chi1inv[i] * D[i];
-          }
-        else
-          for (int iZ = 0; iZ < nZ; iZ++) {
-            const int i = yee_idx + iZ - sR;
-            E[i] = D[i];
-          }
-      }
-    }
-
   return allocated_eh;
 }
 

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;