From 4cb2d363a014eb82db64428bb912be9c9f9a11fd Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 3 Nov 2023 13:59:40 -0700 Subject: [PATCH] =?UTF-8?q?Special=20field=20updates=20for=20r=20=3D=200?= =?UTF-8?q?=20and=20m=20=3D=200,=20=C2=B11=20in=20cylindrical=20coordinate?= =?UTF-8?q?s=20(#2538)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * 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 fc9b1c8fbca5633d55ce685d827a601e278fff16. * 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 --- python/tests/test_adjoint_cyl.py | 2 +- python/tests/test_ldos.py | 17 ++++++-- python/tests/test_pml_cyl.py | 18 +++++---- src/step_db.cpp | 66 ++++++++++++++++++++++++-------- src/update_eh.cpp | 16 +++++++- 5 files changed, 90 insertions(+), 29 deletions(-) diff --git a/python/tests/test_adjoint_cyl.py b/python/tests/test_adjoint_cyl.py index 2995fb614..08a328bb6 100644 --- a/python/tests/test_adjoint_cyl.py +++ b/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) diff --git a/python/tests/test_ldos.py b/python/tests/test_ldos.py index 3b39e4aac..ae77948c1 100644 --- a/python/tests/test_ldos.py +++ b/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, diff --git a/python/tests/test_pml_cyl.py b/python/tests/test_pml_cyl.py index 39ed61c9a..95942f263 100644 --- a/python/tests/test_pml_cyl.py +++ b/python/tests/test_pml_cyl.py @@ -20,7 +20,8 @@ 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), ] @@ -28,14 +29,13 @@ def setUp(cls): 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], diff --git a/src/step_db.cpp b/src/step_db.cpp index 9a54fa5fb..e9a08300d 100644 --- a/src/step_db.cpp +++ b/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]); diff --git a/src/update_eh.cpp b/src/update_eh.cpp index 040806db3..31e1c9e6c 100644 --- a/src/update_eh.cpp +++ b/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]); + } + } + } } } }