Special field updates for r = 0 and m = 0, ±1 in cylindrical coordina…

…tes (#2538)

* special case for update E from D in cyl coords for m = 0 and r = 0

* add new function step_update_EDHB0 in step_generic.cpp which includes PML and nonlinearities

* remove STEP_UPDATE_EDHB0 function and instead call STEP_UPDATE_EDHB with null off-diagonal components

* remove STEP_UPDATE_EDHB0 in step_generic.cpp

* use narrowband source for m=-1 case and add new test for m=0

* add comment to explain why narrow bandwidth source is necessary for |m| = 1

* add step-db update equations for z-PML and absorber at r = 0 for |m| = 1

* Revert "add step-db update equations for z-PML and absorber at r = 0 for |m| = 1"

This reverts commit fc9b1c8.

* Update step_db.cpp

* typos

* fix errors

* remove hack for m = -1 test case

* update r=0 step_db for m=0 to match step_generic

* flip sign of curl terms and strides

* skip special case of r=0 fields check in test_simple_metallic of cylindrical.cpp

* reduce error tolerance in field comparison at r=0 in test_pml of cylindrical.cpp

* check converged flux matches to seven rather than eight digits in test_pml_cyl.py

* revert sign flip in sd of step_db.cpp for |m|=1 update equations

* reduce range of iz-index loop for Dp update at r=0 for |m|=1 by one

* replace custom loops for r=0 update of m = 0, ±1 with LOOP_OVER_IVECS

* compute radiated flux from the farfield radiation pattern rather than Poynting flux of nearfields to avoid DFT interpolation to voxel center

* revert extraction efficiency calculation to using add_flux but with Er source at r > 0

---------

Co-authored-by: Steven G. Johnson <[email protected]>

python/tests/test_adjoint_cyl.py

@@ -202,7 +202,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.3
+        tol = 0.6 if m == 0 else 0.3
         self.assertClose(adj_scale, fd_grad, epsilon=tol)
 
 

python/tests/test_ldos.py

@@ -214,7 +214,13 @@ def ext_eff_cyl(self, dmat, h, m):
         ]
 
         src_cmpt = mp.Er
-        src_pt = mp.Vector3(0, 0, -0.5 * sz + h * dmat)
+
+        # Because (1) Er is not defined at r=0 on the Yee grid, and (2) there
+        # seems to be a bug in the restriction of an Er point source at r = 0,
+        # the source is placed at r = ~Δr (just outside the first voxel).
+        # This incurs a small error which decreases linearly with resolution.
+        src_pt = mp.Vector3(1.5 / self.resolution, 0, -0.5 * sz + h * dmat)
+
         sources = [
             mp.Source(
                 src=mp.GaussianSource(self.fcen, fwidth=0.1 * self.fcen),
@@ -263,7 +269,12 @@ def ext_eff_cyl(self, dmat, h, m):
         )
 
         out_flux = mp.get_fluxes(flux_air)[0]
-        dV = np.pi / (self.resolution**3)
+
+        if src_pt.x == 0:
+            dV = np.pi / (self.resolution**3)
+        else:
+            dV = 2 * np.pi * src_pt.x / (self.resolution**2)
+
         total_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV
         ext_eff = out_flux / total_flux
         print(f"extraction efficiency (cyl):, " f"{dmat:.4f}, {h:.4f}, {ext_eff:.6f}")
@@ -436,7 +447,7 @@ def test_ldos_ext_eff(self):
         ext_eff_cyl = self.ext_eff_cyl(layer_thickness, dipole_height, -1.0)
         ext_eff_3D = self.ext_eff_3D(layer_thickness, dipole_height)
 
-        self.assertAlmostEqual(ext_eff_cyl, ext_eff_3D, delta=0.01)
+        self.assertAlmostEqual(ext_eff_cyl, ext_eff_3D, delta=0.02)
 
         ext_eff_cyl_m_plus = self.ext_eff_cyl(
             layer_thickness,

python/tests/test_pml_cyl.py

@@ -20,22 +20,22 @@ def setUp(cls):
 
     @parameterized.parameterized.expand(
         [
-            (-1.0, 0.0, False),
+            (0.0, 0.04, False),
+            (-1.0, 0, False),
             (2.0, 0.14, False),
             (3.0, 0.17, True),
         ]
     )
     def test_pml_cyl(
         self, m: float, rpos: float, accurate_fields_near_cylorigin: bool = False
     ):
-        """Verifies that the z-PML in cylindrical coordinates properly
-        attenuates fields at r=0.
+        """Verifies that the z-PML properly attenuates fields at r=0.
 
         Args:
            m: exp(imϕ) angular dependence of the fields.
            rpos: position of source along R direction.
            accurate_fields_near_cylorigin: whether to compute more accurate
-              fields near the origin r=0.
+              fields near r=0 for |m| > 1.
         """
 
         pml_layers = [
@@ -119,12 +119,16 @@ def test_pml_cyl(
                 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)
+            cur_flux_str = ", ".join(f"{c:.8f}" for c in cur_flux)
             flux_tot = sum(cur_flux)
 
-            print(f"flux:, {sim.meep_time()}, {cur_flux_str}, {flux_tot:.6f}")
+            print(f"flux:, {sim.meep_time()}, {cur_flux_str}, {flux_tot:.8f}")
 
-            places = 6 if mp.is_single_precision() else 8
+            # Check that the flux is converged with runtime long after the
+            # source has turned off. This verifies the correctness of the
+            # z-PML at r=0 for m=0, ±1 which involve special field-update
+            # equations.
+            places = 6 if mp.is_single_precision() else 7
             for i in range(len(cur_flux)):
                 self.assertAlmostEqual(
                     prev_flux[i],

src/step_db.cpp

@@ -281,21 +281,38 @@ bool fields_chunk::step_db(field_type ft) {
         // d(Dz)/dt = (1/r) * d(r*Hp)/dr
         const realnum *g = f[Hp][cmp];
         const realnum *cndinv = s->condinv[Dz][Z];
+        const realnum *cnd = s->conductivity[Dz][Z];
         realnum *fcnd = f_cond[Dz][cmp];
         const direction dsig = cycle_direction(gv.dim, Z, 1);
         const realnum *siginv = s->sigsize[dsig] > 1 ? s->siginv[dsig] : 0;
-        const int dk = gv.iyee_shift(Dz).in_direction(dsig);
+        const realnum *sig = s->sigsize[dsig] > 1 ? s->sig[dsig] : 0;
+        const realnum *kap = s->sigsize[dsig] > 1 ? s->kap[dsig] : 0;
         const direction dsigu = cycle_direction(gv.dim, Z, 2);
         const realnum *siginvu = s->sigsize[dsigu] > 1 ? s->siginv[dsigu] : 0;
-        const int dku = gv.iyee_shift(Dz).in_direction(dsigu);
+        const realnum *sigu = s->sigsize[dsigu] > 1 ? s->sig[dsigu] : 0;
+        const realnum *kapu = s->sigsize[dsigu] > 1 ? s->kap[dsigu] : 0;
         realnum *fu = siginvu && f_u[Dz][cmp] ? f[Dz][cmp] : 0;
         realnum *the_f = fu ? f_u[Dz][cmp] : f[Dz][cmp];
-        for (int iz = 0; iz < nz; ++iz) {
-          // Note: old code (prior to Meep 0.2) was missing factor of 4??
-          realnum df, dfcnd = g[iz] * (Courant * 4) * (cndinv ? cndinv[iz] : 1);
-          if (fcnd) fcnd[iz] += dfcnd;
-          the_f[iz] += (df = dfcnd * (siginv ? siginv[dk + 2 * (dsig == Z) * iz] : 1));
-          if (fu) fu[iz] += siginvu[dku + 2 * (dsigu == Z) * iz] * df;
+        realnum dt2 = dt * 0.5;
+
+        ivec is = gv.little_owned_corner(Dz);
+        ivec ie = gv.big_owned_corner(Dz);
+        ie.set_direction(R, 0);
+        LOOP_OVER_IVECS(gv, is, ie, i) {
+          realnum fprev = the_f[i];
+          realnum dfcnd = g[i] * (Courant * 4);
+          if (fcnd) {
+            realnum fcnd_prev = fcnd[i];
+            fcnd[i] = ((1 - dt2 * cnd[i]) * fcnd[i] + dfcnd) * cndinv[i];
+            dfcnd = fcnd[i] - fcnd_prev;
+          }
+          KSTRIDE_DEF(dsig, k, is, gv);
+          DEF_k;
+          KSTRIDE_DEF(dsigu, ku, is, gv);
+          DEF_ku;
+          the_f[i] = ((kap ? kap[k] - sig[k] : 1) * the_f[i] + dfcnd) * (siginv ? siginv[k] : 1);
+          if (fu)
+            fu[i] = siginvu[ku] * ((kapu ? kapu[ku] - sigu[ku] : 1) * fu[i] + the_f[i] - fprev);
         }
         ZERO_Z(f[Dp][cmp]);
         if (f_cond[Dp][cmp]) ZERO_Z(f_cond[Dp][cmp]);
@@ -315,25 +332,40 @@ bool fields_chunk::step_db(field_type ft) {
         const realnum *f_p = f[ft == D_stuff ? Hr : Ep][cmp];
         const realnum *f_m = ft == D_stuff ? f[Hz][cmp] : (f[Ez][1 - cmp] + (nz + 1));
         const realnum *cndinv = s->condinv[cc][d_c];
+        const realnum *cnd = s->conductivity[cc][d_c];
         realnum *fcnd = f_cond[cc][cmp];
         const direction dsig = cycle_direction(gv.dim, d_c, 1);
         const realnum *siginv = s->sigsize[dsig] > 1 ? s->siginv[dsig] : 0;
-        const int dk = gv.iyee_shift(cc).in_direction(dsig);
+        const realnum *sig = s->sigsize[dsig] > 1 ? s->sig[dsig] : 0;
+        const realnum *kap = s->sigsize[dsig] > 1 ? s->kap[dsig] : 0;
         const direction dsigu = cycle_direction(gv.dim, d_c, 2);
         const realnum *siginvu = s->sigsize[dsigu] > 1 ? s->siginv[dsigu] : 0;
-        const int dku = gv.iyee_shift(cc).in_direction(dsigu);
+        const realnum *sigu = s->sigsize[dsigu] > 1 ? s->sig[dsigu] : 0;
+        const realnum *kapu = s->sigsize[dsigu] > 1 ? s->kap[dsigu] : 0;
         realnum *fu = siginvu && f_u[cc][cmp] ? f[cc][cmp] : 0;
         realnum *the_f = fu ? f_u[cc][cmp] : f[cc][cmp];
         int sd = ft == D_stuff ? +1 : -1;
         realnum f_m_mult = ft == D_stuff ? 2 : (1 - 2 * cmp) * m;
+        realnum dt2 = dt * 0.5;
 
-        for (int iz = (ft == D_stuff); iz < nz + (ft == D_stuff); ++iz) {
-          realnum df;
-          realnum dfcnd = (sd * Courant) * (f_p[iz] - f_p[iz - sd] - f_m_mult * f_m[iz]) *
-                          (cndinv ? cndinv[iz] : 1);
-          if (fcnd) fcnd[iz] += dfcnd;
-          the_f[iz] += (df = dfcnd * (siginv ? siginv[dk + 2 * (dsig == Z) * iz] : 1));
-          if (fu) fu[iz] += siginvu[dku + 2 * (dsigu == Z) * iz] * df;
+        ivec is = gv.little_owned_corner(cc);
+        ivec ie = gv.big_owned_corner(cc);
+        ie.set_direction(R, 0);
+        LOOP_OVER_IVECS(gv, is, ie, i) {
+          realnum fprev = the_f[i];
+          realnum dfcnd = (sd * Courant) * (f_p[i] - f_p[i - sd] - f_m_mult * f_m[i]);
+          if (fcnd) {
+            realnum fcnd_prev = fcnd[i];
+            fcnd[i] = ((1 - dt2 * cnd[i]) * fcnd[i] + dfcnd) * cndinv[i];
+            dfcnd = fcnd[i] - fcnd_prev;
+          }
+          KSTRIDE_DEF(dsig, k, is, gv);
+          DEF_k;
+          KSTRIDE_DEF(dsigu, ku, is, gv);
+          DEF_ku;
+          the_f[i] = ((kap ? kap[k] - sig[k] : 1) * the_f[i] + dfcnd) * (siginv ? siginv[k] : 1);
+          if (fu)
+            fu[i] = siginvu[ku] * ((kapu ? kapu[ku] - sigu[ku] : 1) * fu[i] + the_f[i] - fprev);
         }
         if (ft == D_stuff) {
           ZERO_Z(f[Dz][cmp]);

src/update_eh.cpp

@@ -187,12 +187,26 @@ bool fields_chunk::update_eh(field_type ft, bool skip_w_components) {
                  sizeof(realnum) * gv.ntot());
         }
 
-        if (f[ec][cmp] != f[dc][cmp])
+        if (f[ec][cmp] != f[dc][cmp]) {
           STEP_UPDATE_EDHB(f[ec][cmp], ec, gv, gvs_eh[ft][i].little_owned_corner0(ec),
                            gvs_eh[ft][i].big_corner(), dmp[dc][cmp], dmp[dc_1][cmp], dmp[dc_2][cmp],
                            s->chi1inv[ec][d_ec], dmp[dc_1][cmp] ? s->chi1inv[ec][d_1] : NULL,
                            dmp[dc_2][cmp] ? s->chi1inv[ec][d_2] : NULL, s_ec, s_1, s_2, s->chi2[ec],
                            s->chi3[ec], f_w[ec][cmp], dsigw, s->sig[dsigw], s->kap[dsigw]);
+
+          if (gv.dim == Dcyl) {
+            ivec is = gvs_eh[ft][i].little_owned_corner(ec);
+            if (is.r() == 0) {
+              ivec ie = gvs_eh[ft][i].big_corner();
+              ie.set_direction(R, 0);
+              /* pass NULL for off-diagonal terms since they must be
+                 zero at r=0 for an axisymmetric structure: */
+              STEP_UPDATE_EDHB(f[ec][cmp], ec, gv, is, ie, dmp[dc][cmp], NULL, NULL,
+                               s->chi1inv[ec][d_ec], NULL, NULL, s_ec, s_1, s_2, s->chi2[ec],
+                               s->chi3[ec], f_w[ec][cmp], dsigw, s->sig[dsigw], s->kap[dsigw]);
+            }
+          }
+        }
       }
     }
   }