From 157d37ead3804e6ef4f5956ccb7f187c36a23fe5 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Mon, 31 Aug 2026 11:18:32 -0700 Subject: [PATCH 1/3] Route sigma through a volume-aware hook `structure_chunk::set_chi1inv` volume-averages the instantaneous permittivity while `add_susceptibility` point-samples sigma, so a dispersive material switches on a whole pixel at a time as an object moves. That is a staircase, and a staircase has no derivative -- the same wall the geometry gradients hit, reached through the material rather than the shape. Add `material_function::eff_sigma_row`, standing to `sigma_row` as `eff_chi1inv_row` stands to `chi1p1`, and call it from the two sigma loops with the pixel volume the chi1inv loop already uses. The default implementation samples the volume centre, so this is behaviour-preserving on its own: nothing overrides it yet. The sigma loops now index in ivec like the chi1inv loop above them, since the half-pixel offdiagonal shift and the smoothing diameter have to agree between the two for the eventual averaging to be consistent. Verified neutral: a Drude MaterialGrid's transmittance is unchanged bit for bit, and test_adjoint_dispersion, test_material_grid, test_material_dispersion, test_dispersive_eigenmode, test_faraday_rotation, test_multilevel_atom, test_geometry_gradient and test_conductivity all pass. --- src/anisotropic_averaging.cpp | 17 ++++++++++------- src/meep.hpp | 12 ++++++++++++ 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/src/anisotropic_averaging.cpp b/src/anisotropic_averaging.cpp index 7ce588510..bef4387f0 100644 --- a/src/anisotropic_averaging.cpp +++ b/src/anisotropic_averaging.cpp @@ -361,7 +361,10 @@ void structure_chunk::add_susceptibility(material_function &sigma, field_type ft realnum *s0 = newsus->sigma[c][d0]; realnum *s1 = newsus->sigma[c][d1]; realnum *s2 = newsus->sigma[c][d2]; - vec shift1(gv[unit_ivec(gv.dim, component_direction(c)) * (ft == E_stuff ? 1 : -1)]); + // Same half-pixel offdiagonal convention, and same one-pixel smoothing + // diameter, as the chi1inv loop in set_chi1inv() above. + const double smoothing_diameter = 1.0; + ivec shift1(unit_ivec(gv.dim, component_direction(c)) * (ft == E_stuff ? 1 : -1)); // Use OpenMP parallelization when the material function is thread-safe // (i.e., pure C++ geometry, not a Python callback). The trivial[] flags // need a reduction since each thread computes its local subset. @@ -371,9 +374,9 @@ void structure_chunk::add_susceptibility(material_function &sigma, field_type ft gv.big_corner() + gv.iyee_shift(c), i, "omp parallel for collapse(3) reduction(&:trivial0,trivial1,trivial2)") { double sigrow[3], sigrow_offdiag[3]; - IVEC_LOOP_LOC(gv, here); - sigma.sigma_row(c, sigrow, here); - sigma.sigma_row(c, sigrow_offdiag, here - shift1); + IVEC_LOOP_ILOC(gv, here); + sigma.eff_sigma_row(c, sigrow, gv.dV(here, smoothing_diameter)); + sigma.eff_sigma_row(c, sigrow_offdiag, gv.dV(here - shift1, smoothing_diameter)); sigrow[(idiag + 1) % 3] = sigrow_offdiag[(idiag + 1) % 3]; sigrow[(idiag + 2) % 3] = sigrow_offdiag[(idiag + 2) % 3]; if (s0) trivial0 &= (s0[i] = sigrow[0]) == 0.; @@ -388,9 +391,9 @@ void structure_chunk::add_susceptibility(material_function &sigma, field_type ft // Serial path for non-thread-safe material functions (Python callbacks) LOOP_OVER_VOL(gv, c, i) { double sigrow[3], sigrow_offdiag[3]; - IVEC_LOOP_LOC(gv, here); - sigma.sigma_row(c, sigrow, here); - sigma.sigma_row(c, sigrow_offdiag, here - shift1); + IVEC_LOOP_ILOC(gv, here); + sigma.eff_sigma_row(c, sigrow, gv.dV(here, smoothing_diameter)); + sigma.eff_sigma_row(c, sigrow_offdiag, gv.dV(here - shift1, smoothing_diameter)); sigrow[(idiag + 1) % 3] = sigrow_offdiag[(idiag + 1) % 3]; sigrow[(idiag + 2) % 3] = sigrow_offdiag[(idiag + 2) % 3]; if (s0 && (s0[i] = sigrow[0]) != 0.) trivial[0] = false; diff --git a/src/meep.hpp b/src/meep.hpp index 94f0c2c4e..d67d5f1d3 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -550,6 +550,18 @@ class material_function { sigrow[0] = sigrow[1] = sigrow[2] = 0.0; } + /* Volume-averaged sigma over a pixel, standing to sigma_row as + eff_chi1inv_row stands to chi1p1. The instantaneous permittivity has been + subpixel-averaged since forever while sigma was point sampled, which makes + a dispersive material switch on a whole pixel at a time as an object moves + -- a staircase, and a staircase has no derivative. + + The default reproduces the point sample exactly, so a subclass that does + not override this behaves as it always has. */ + virtual void eff_sigma_row(component c, double sigrow[3], const volume &v) { + sigma_row(c, sigrow, v.center()); + } + // Nonlinear susceptibilities virtual bool has_chi3(component c) { (void)c; From d87ec20b515bcc2f7e3fc5ca0a7939bd1d52d1d7 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Mon, 31 Aug 2026 12:15:13 -0700 Subject: [PATCH 2/3] Subpixel-average sigma across a geometric interface Taper each medium's pole amplitude by the pixel's filling fraction, using the same front object, shift and fill that eff_chi1inv_matrix uses, so the permittivity and the susceptibility describe the same interface. Swept sub-pixel across one pixel, a Drude slab whose eps_infinity matches the background -- so that the only position dependence there is comes from sigma: point sampled constant to the last digit at 10 of 11 offsets, then a single 20% jump averaged smooth, symmetric about the pixel centre, and 60x less total variation across the sweep Linear in fill rather than Kottke, deliberately. Kottke's normal component is a harmonic average of a permittivity; there is no analogous exact result for a resonant amplitude, and a harmonic average against a medium that lacks the pole would be zero in every boundary pixel. Scaling each medium's whole sigma tensor by a non-negative weight and adding keeps the result positive semidefinite, which is what passivity depends on, and makes the amplitude exactly linear in fill, so the shape derivative is a constant rather than something to be differenced. Multilevel atoms and gyrotropic media are excluded and stay point sampled. A multilevel atom carries per-pixel level populations with their own nonlinear dynamics, so a partially filled pixel is not a weaker gain medium -- without this exclusion test_multilevel_atom fails. A gyrotropic sigma is antisymmetric, so the positive-semidefiniteness argument does not hold for it either. Material grids and user materials also stay point sampled: their material already varies continuously in space, and mixing a level set with a geometric filling fraction is a separate question. The pole lookup is factored out of sigma_row into material_sigma_row so both paths ask the same question of a material. --- python/Makefile.am | 1 + python/tests/test_dispersive_smoothing.py | 209 ++++++++++++++++++++++ src/meepgeom.cpp | 142 +++++++++++---- src/meepgeom.hpp | 2 + 4 files changed, 316 insertions(+), 38 deletions(-) create mode 100644 python/tests/test_dispersive_smoothing.py diff --git a/python/Makefile.am b/python/Makefile.am index a72c154ac..59994d26f 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -64,6 +64,7 @@ TESTS = \ $(TEST_DIR)/test_cyl_ellipsoid.py \ $(TEST_DIR)/test_dft_energy.py \ $(TEST_DIR)/test_dft_fields.py \ + $(TEST_DIR)/test_dispersive_smoothing.py \ $(DIFFRACTED_PLANEWAVE_TEST) \ $(DISPERSIVE_EIGENMODE_TEST) \ $(TEST_DIR)/test_divide_mpi_processes.py \ diff --git a/python/tests/test_dispersive_smoothing.py b/python/tests/test_dispersive_smoothing.py new file mode 100644 index 000000000..4defbcb60 --- /dev/null +++ b/python/tests/test_dispersive_smoothing.py @@ -0,0 +1,209 @@ +"""Subpixel averaging of susceptibility amplitudes at a geometric interface. + +The instantaneous permittivity has been subpixel-smoothed since forever while +sigma was point sampled, so a dispersive object's response switched on a whole +pixel at a time as it moved. `eff_sigma_row` tapers each medium's pole +amplitude by the pixel's filling fraction instead. + +Every medium here has `eps_infinity = 1`, matching the background, so the +instantaneous permittivity is uniform and contributes no position dependence at +all: the whole signal in these tests comes from sigma. `eps_averaging=False` +routes back to the point sample, so the before/after comparison comes out of +one binary. +""" + +import unittest + +import numpy as np + +import meep as mp + +RES = 20 +FCEN = 1.0 +RUN = 120.0 +DX = 1.0 / RES + +DRUDE = mp.Medium( + epsilon=1.0, + E_susceptibilities=[mp.DrudeSusceptibility(frequency=1.0, gamma=0.5, sigma=1.0)], +) +LOSSLESS = mp.Medium( + epsilon=1.0, + E_susceptibilities=[ + mp.LorentzianSusceptibility(frequency=1.6, gamma=0.0, sigma=0.6) + ], +) + + +def transmission(offset, averaging=True, res=RES, material=DRUDE, background=mp.air): + sim = mp.Simulation( + cell_size=mp.Vector3(6, 4), + resolution=res, + boundary_layers=[mp.PML(1.0)], + default_material=background, + sources=[ + mp.Source( + mp.GaussianSource(FCEN, fwidth=0.2), + component=mp.Ez, + center=mp.Vector3(-1.8, 0), + ) + ], + geometry=[ + mp.Block( + center=mp.Vector3(offset, 0), + size=mp.Vector3(0.5, mp.inf), + material=material, + ) + ], + eps_averaging=averaging, + force_complex_fields=True, + ) + mon = sim.add_dft_fields( + [mp.Ez], FCEN, 0, 1, center=mp.Vector3(1.8, 0), size=mp.Vector3() + ) + sim.run(until=RUN) + return float(np.abs(sim.get_dft_array(mon, mp.Ez, 0)) ** 2) + + +def step_ratio(curve): + """Largest single step relative to the average step. + + A staircase puts all of its variation in one place, so this is close to the + number of samples; a smooth curve sampled this finely keeps it near 1. + """ + steps = np.abs(np.diff(np.asarray(curve))) + return float(steps.max() / steps.mean()) if steps.mean() > 0 else np.inf + + +class TestSubpixelContinuity(unittest.TestCase): + OFFSETS = np.linspace(0, DX, 11) + + def setUp(self): + mp.verbosity(0) + + def test_point_sampled_sigma_is_a_staircase(self): + # Not a claim about desired behaviour -- it pins what the comparison + # below is measuring against, so a change that made both paths smooth + # for some unrelated reason could not pass vacuously. + curve = [transmission(o, averaging=False) for o in self.OFFSETS] + identical = sum(a == b for a, b in zip(curve, curve[1:])) + self.assertGreaterEqual( + identical, len(curve) - 3, f"expected a staircase, got {curve}" + ) + + def test_averaged_sigma_varies_continuously(self): + curve = [transmission(o, averaging=True) for o in self.OFFSETS] + self.assertEqual( + sum(a == b for a, b in zip(curve, curve[1:])), + 0, + f"adjacent offsets should differ: {curve}", + ) + # Measured ~1.9 for the averaged curve and ~10 (i.e. one step carries + # everything) for the point-sampled one. + self.assertLess(step_ratio(curve), 3.0, f"curve looks stepped: {curve}") + + def test_averaging_shrinks_the_discretization_sensitivity(self): + # The structure is the same slab throughout; all variation across the + # sweep is discretization error. Smoothing should reduce it a lot. + span = lambda c: max(c) - min(c) + rough = span([transmission(o, averaging=False) for o in self.OFFSETS]) + smooth = span([transmission(o, averaging=True) for o in self.OFFSETS]) + self.assertLess(smooth, rough / 10, f"span {smooth} vs {rough}") + + def test_translating_by_a_whole_pixel_reproduces_the_structure(self): + # Exact invariant: shifting by one pixel maps the grid onto itself, so + # any discrepancy is a bug in the fill rather than discretization. + a, b = transmission(0.0), transmission(DX) + self.assertAlmostEqual(a / b, 1.0, places=5) + + +class TestExclusions(unittest.TestCase): + """Pole kinds whose amplitude a filling fraction cannot meaningfully scale. + + A multilevel atom carries per-pixel level populations with their own + nonlinear dynamics, so a partially filled pixel is not a weaker gain + medium. A gyrotropic sigma is antisymmetric, so the argument that a + non-negative combination of two tensors stays positive semidefinite -- the + reason the linear mixture is safe -- does not apply. Both stay point + sampled, and these pin that rather than leaving it to a comment. + """ + + def setUp(self): + mp.verbosity(0) + + def assertStillStepped(self, susceptibility): + medium = mp.Medium(epsilon=1.0, E_susceptibilities=[susceptibility]) + curve = [ + transmission(o, averaging=True, material=medium) + for o in (0.1 * DX, 0.3 * DX, 0.5 * DX) + ] + self.assertEqual( + curve[0], curve[1], f"pole should not have been smoothed: {curve}" + ) + self.assertEqual(curve[1], curve[2], f"pole should not been smoothed: {curve}") + + def test_multilevel_atom_is_not_smoothed(self): + self.assertStillStepped( + mp.MultilevelAtom( + sigma=1, + transitions=[ + mp.Transition(1, 2, pumping_rate=0.005, frequency=FCEN, gamma=1e-5) + ], + initial_populations=[1], + ) + ) + + def test_gyrotropic_is_not_smoothed(self): + self.assertStillStepped( + mp.GyrotropicLorentzianSusceptibility( + frequency=1.1, gamma=0.2, sigma=0.5, bias=mp.Vector3(0, 0, 1) + ) + ) + + +class TestInvariants(unittest.TestCase): + def setUp(self): + mp.verbosity(0) + + def test_a_block_of_the_background_medium_is_invisible(self): + # Exact, and independent of any tolerance: an object made of what is + # already there cannot change the answer no matter where it sits. + a = transmission(0.0, material=DRUDE, background=DRUDE) + b = transmission(0.4 * DX, material=DRUDE, background=DRUDE) + self.assertAlmostEqual(a / b, 1.0, places=9) + + def test_a_lossless_inclusion_does_not_pump_energy(self): + # Passivity, end to end. Mixing two media's sigma tensors is only safe + # if the result stays positive semidefinite; scaling each by a + # non-negative fill and adding does, but this is the check that would + # catch it if the construction were ever changed to something that + # does not. A closed cavity with gamma = 0 has nowhere to lose energy, + # so the field must not grow without bound. + sim = mp.Simulation( + cell_size=mp.Vector3(3, 3), + resolution=RES, + geometry=[ + mp.Block( + center=mp.Vector3(0.4 * DX, 0.4 * DX), + size=mp.Vector3(1.0, 1.0), + material=LOSSLESS, + ) + ], + sources=[ + mp.Source( + mp.GaussianSource(FCEN, fwidth=0.4), + component=mp.Ez, + center=mp.Vector3(-0.7, 0.3), + ) + ], + eps_averaging=True, + ) + sim.run(until=50) + early = sim.field_energy_in_box(mp.Volume(mp.Vector3(), size=mp.Vector3(3, 3))) + sim.run(until=400) + late = sim.field_energy_in_box(mp.Volume(mp.Vector3(), size=mp.Vector3(3, 3))) + self.assertLess(late, 1.05 * early, f"energy grew: {early} -> {late}") + + +if __name__ == "__main__": + unittest.main() diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index b2ad4052a..b0af48cb1 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1698,6 +1698,109 @@ static bool susceptibility_equiv(const susceptibility &o0, const susceptibility return true; } +/* Sum the amplitudes of `mat`'s poles equivalent to the one being registered, + into the c'th row of its sigma tensor. Factored out of sigma_row() so that + eff_sigma_row() can ask the same question of two materials at once. */ +void geom_epsilon::material_sigma_row(meep::component c, double sigrow[3], material_type mat) { + sigrow[0] = sigrow[1] = sigrow[2] = 0.0; + if (!current_pol) return; + + if (mat->which_subclass != material_data::MATERIAL_USER && + mat->which_subclass != material_data::MATERIAL_GRID && + mat->which_subclass != material_data::MEDIUM) + return; + + const susceptibility_list &slist = + type(c) == meep::E_stuff ? mat->medium.E_susceptibilities : mat->medium.H_susceptibilities; + /* Accumulate rather than stopping at the first match. A material grid's + list is the concatenation of its two constituents' poles, so when both + media carry the same pole -- same frequency and damping, different + strength -- two entries are equivalent to the one being registered and + the mixture's amplitude is their sum. This is also the right answer for + a plain medium that happens to list a pole twice, since identical poles + superpose. */ + for (const susceptibility &susc : slist) { + if (!susceptibility_equiv(susc, current_pol->user_s)) continue; + switch (meep::component_index(c)) { // which row of the sigma tensor to return + case 0: + sigrow[0] += susc.sigma_diag.x; + sigrow[1] += susc.sigma_offdiag.x; + sigrow[2] += susc.sigma_offdiag.y; + break; + case 1: + sigrow[0] += susc.sigma_offdiag.x; + sigrow[1] += susc.sigma_diag.y; + sigrow[2] += susc.sigma_offdiag.z; + break; + default: // case 2: + sigrow[0] += susc.sigma_offdiag.y; + sigrow[1] += susc.sigma_offdiag.z; + sigrow[2] += susc.sigma_diag.z; + break; + } + } +} + +/* Volume-average sigma across a geometric interface. + + The permittivity has always been subpixel-smoothed while sigma was point + sampled, so a dispersive object's response switched on a whole pixel at a + time as it moved. Taper each medium's pole amplitude by the pixel's filling + fraction instead, using the same front object, shift and fill that + eff_chi1inv_matrix uses, so the two describe the same interface. + + Linear in fill rather than Kottke, deliberately. Kottke's normal component + is a harmonic average of a permittivity; there is no analogous exact result + for a resonant amplitude, and a harmonic average against a medium without + the pole would be zero in every boundary pixel. Scaling each medium's whole + sigma tensor by a non-negative weight and adding keeps the result positive + semidefinite -- which is what passivity depends on -- and makes the + amplitude exactly linear in fill, so the shape derivative is a constant. + + Material grids and user materials are point sampled as before: their + material varies continuously in space already, and mixing a level set with a + geometric filling fraction is a separate question. */ +void geom_epsilon::eff_sigma_row(meep::component c, double sigrow[3], const meep::volume &v) { + const geometric_object *o; + material_type mat, mat_behind; + vector3 p, shiftby; + + // No smoothing requested: match what the permittivity is doing. + if (maxeval == 0 || !current_pol) return sigma_row(c, sigrow, v.center()); + + /* Out of scope, point sample instead. A multilevel atom carries per-pixel + level populations with their own nonlinear dynamics, so there is no reason + to think scaling an amplitude by a filling fraction gives the response of + a partially filled pixel the way it does for a linear pole. A gyrotropic + sigma is antisymmetric, so the argument that a non-negative combination + stays positive semidefinite does not apply to it either. */ + const susceptibility &pole = current_pol->user_s; + if (!pole.transitions.empty() || !pole.initial_populations.empty() || pole.saturated_gyrotropy || + pole.bias.x != 0 || pole.bias.y != 0 || pole.bias.z != 0) + return sigma_row(c, sigrow, v.center()); + + if (!get_front_object(v, geometry_tree, p, &o, shiftby, mat, mat_behind)) + return sigma_row(c, sigrow, v.center()); + + if (material_type_equal(mat, mat_behind) || is_variable(mat) || is_variable(mat_behind)) { + material_gc(mat); + return sigma_row(c, sigrow, v.center()); + } + + geom_box pixel = gv2box(v); + pixel.low = vector3_minus(pixel.low, shiftby); + pixel.high = vector3_minus(pixel.high, shiftby); + double fill = box_overlap_with_object(pixel, *o, tol, maxeval); + + double front[3], behind[3]; + material_sigma_row(c, front, mat); + material_sigma_row(c, behind, mat_behind); + for (int i = 0; i < 3; ++i) + sigrow[i] = fill * front[i] + (1 - fill) * behind[i]; + + material_gc(mat); +} + void geom_epsilon::sigma_row(meep::component c, double sigrow[3], const meep::vec &r) { vector3 p = vec_to_vector3(r); @@ -1724,44 +1827,7 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3], const meep::ve mat->medium.check_offdiag_im_zero_or_abort(); } - sigrow[0] = sigrow[1] = sigrow[2] = 0.0; - - if (mat->which_subclass == material_data::MATERIAL_USER || - mat->which_subclass == material_data::MATERIAL_GRID || - mat->which_subclass == material_data::MEDIUM) { - - const susceptibility_list &slist = - type(c) == meep::E_stuff ? mat->medium.E_susceptibilities : mat->medium.H_susceptibilities; - /* Accumulate rather than stopping at the first match. A material grid's - list is the concatenation of its two constituents' poles, so when both - media carry the same pole -- same frequency and damping, different - strength -- two entries are equivalent to the one being registered and - the mixture's amplitude is their sum. This is also the right answer for - a plain medium that happens to list a pole twice, since identical poles - superpose. */ - for (const susceptibility &susc : slist) { - if (susceptibility_equiv(susc, current_pol->user_s)) { - int ic = meep::component_index(c); - switch (ic) { // which row of the sigma tensor to return - case 0: - sigrow[0] += susc.sigma_diag.x; - sigrow[1] += susc.sigma_offdiag.x; - sigrow[2] += susc.sigma_offdiag.y; - break; - case 1: - sigrow[0] += susc.sigma_offdiag.x; - sigrow[1] += susc.sigma_diag.y; - sigrow[2] += susc.sigma_offdiag.z; - break; - default: // case 2: - sigrow[0] += susc.sigma_offdiag.y; - sigrow[1] += susc.sigma_offdiag.z; - sigrow[2] += susc.sigma_diag.z; - break; - } - } - } - } + material_sigma_row(c, sigrow, mat); material_gc(mat); } diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index d8f0c6777..5f485e9cd 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -239,6 +239,8 @@ class geom_epsilon : public meep::material_function { double tol, int maxeval); virtual void sigma_row(meep::component c, double sigrow[3], const meep::vec &r); + virtual void eff_sigma_row(meep::component c, double sigrow[3], const meep::volume &v); + void material_sigma_row(meep::component c, double sigrow[3], material_type mat); void add_susceptibilities(meep::structure *s); void add_susceptibilities(meep::field_type ft, meep::structure *s); From 373fc6031400922ba7c8a4138b8743a6d0dfc9f2 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Mon, 31 Aug 2026 13:08:46 -0700 Subject: [PATCH 3/3] Shape derivatives for dispersive objects `geometry_addgradient` contracted only the real, non-dispersive `eff_chi1inv_matrix`, so a Drude or Lorentz object took a branch that dropped its poles entirely: the gradient described a structure the forward solve never ran. That is the piece a metal reflector's position needs. Add `eff_chi1inv_row_disp_fill`, the dispersive counterpart of the `fill_override` on `eff_chi1inv_matrix`. `eff_chi1inv_row_disp` cannot serve, because it resolves the material at a point and so knows nothing about a filling fraction. The construction matches what the forward solve now runs: the instantaneous part is Kottke-smoothed, and each pole's amplitude is taken linearly in fill from whichever medium carries it, exactly as `eff_sigma_row` writes into the sigma arrays. Since chi1 is linear in sigma, a pole contributes `lineshape(freq) * weight * sigma` and the two media's pole lists never have to be matched against each other -- the same reason `epsilon_material_grid` can concatenate rather than blend. The branch is chosen from the materials at *this* interface rather than from `md->trivial` of the front object. That distinction matters: `trivial` is a property of a material, but under mixing it becomes a property of a pixel, so a plain dielectric sitting in front of a metal would otherwise take the non-dispersive branch and drop the whole contribution at precisely the boundary pixels that carry the derivative. `interface_fill` now hands back both materials, since the caller needs them to build the mixture. Against a central difference of the objective, at resolution 20: dielectric (control) 0.139% Drude metal 0.017% both inside the 2% the existing shape-derivative tests use, which is set by the O(dx) accuracy of a shape derivative rather than by anything dispersive. --- python/tests/test_geometry_gradient.py | 42 +++++++ src/meepgeom.cpp | 155 +++++++++++++++++++++---- src/meepgeom.hpp | 3 +- 3 files changed, 175 insertions(+), 25 deletions(-) diff --git a/python/tests/test_geometry_gradient.py b/python/tests/test_geometry_gradient.py index 3a3d24595..f20d10412 100644 --- a/python/tests/test_geometry_gradient.py +++ b/python/tests/test_geometry_gradient.py @@ -254,6 +254,48 @@ def value(dy): self.assertGreater(abs(reference), 1e-6, "objective must respond to the move") self.assertLess(abs(adjoint - reference) / abs(reference), 2e-2) + def test_dispersive_object(self): + # A metal reflector's position is the motivating case. Before the + # dispersive branch was routed, `geometry_addgradient` contracted only + # the real non-dispersive tensor, so a Drude block's poles contributed + # nothing and the gradient was of a structure that was not simulated. + drude = mp.Medium( + epsilon=1.0, + E_susceptibilities=[ + mp.DrudeSusceptibility(frequency=1.0, gamma=0.2, sigma=0.6) + ], + ) + size = mp.Vector3(1.0, 0.3) + offset = quarter_pixel(RES) + + def block(center, differentiable=None): + return mp.Block( + center=mp.Vector3(center.x, center.y + offset), + size=mp.Vector3(size.x + 2 * offset, size.y), + material=drude, + differentiable=differentiable, + name="scatterer", + ) + + opt = _problem(block(mp.Vector3(), differentiable=["center"])) + _, grad = opt([RHO]) + adjoint = float(np.real(np.atleast_1d(grad["scatterer"]["center"])[1])) + + def value(dy): + return float( + np.asarray( + _problem(block(mp.Vector3(0, dy)))([RHO], need_gradient=False)[0] + ).real.item() + ) + + reference = (value(FD_STEP) - value(-FD_STEP)) / (2 * FD_STEP) + self.assertGreater(abs(reference), 1e-6, "objective must respond to the move") + self.assertLess( + abs(adjoint - reference) / abs(reference), + 2e-2, + f"adjoint {adjoint} vs finite difference {reference}", + ) + def test_no_flagged_object_leaves_the_return_shape_alone(self): opt = _problem(_block(mp.Vector3(), self.SIZE)) _, gradient = opt([RHO]) diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index b0af48cb1..be776eeed 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1221,7 +1221,8 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_ma which case the shape derivative vanishes there and the caller can skip it -- which is what restricts the gradient loop to the boundary shell. */ bool geom_epsilon::interface_fill(const meep::volume &v, double tol, int maxeval, double &fill, - const geometric_object **which, vector3 *shift) { + const geometric_object **which, vector3 *shift, + material_type *front_mat, material_type *behind_mat) { const geometric_object *o; material_type mat, mat_behind; vector3 p, shiftby; @@ -1236,6 +1237,8 @@ bool geom_epsilon::interface_fill(const meep::volume &v, double tol, int maxeval fill = box_overlap_with_object(pixel, *o, tol, maxeval); if (which) *which = o; if (shift) *shift = shiftby; + if (front_mat) *front_mat = mat; + if (behind_mat) *behind_mat = mat_behind; return true; } @@ -2745,6 +2748,86 @@ void get_chi1_tensor_disp(std::complex tensor[9], const meep::vec &r, do } } +/* Does this material contribute anything the non-dispersive tensor misses? */ +bool is_dispersive(material_type m) { + if (!m) return false; + const medium_struct *mm = &m->medium; + return !mm->E_susceptibilities.empty() || mm->D_conductivity_diag.x != 0 || + mm->D_conductivity_diag.y != 0 || mm->D_conductivity_diag.z != 0; +} + +/* The effective dispersive tensor of a pixel straddling an interface, at a + filling fraction supplied by the caller. + + This is the dispersive counterpart of the `fill_override` on + eff_chi1inv_matrix, and it exists for the same reason: a shape derivative + needs d(chi1inv)/d(fill), and getting it from meep's own assembly is better + than re-deriving it. eff_chi1inv_row_disp cannot serve, because it resolves + the material at a *point* and so knows nothing about a filling fraction. + + The construction matches what the forward solve is actually running. The + instantaneous part is Kottke-smoothed, exactly as structure_chunk::set_chi1inv + does it; each pole's amplitude is taken linearly in fill from whichever + medium carries it, exactly as eff_sigma_row now writes into the sigma arrays. + Because chi1 is linear in sigma, a pole contributes L(freq) * weight * sigma + and the two media's pole lists never have to be matched against each other -- + the same reason epsilon_material_grid can concatenate rather than blend. */ +static bool eff_chi1inv_row_disp_fill(meep::component c, std::complex chi1inv_row[3], + const meep::volume &v, double freq, geom_epsilon *geps, + double fill, material_type front, material_type behind) { + bool fallback = false; + symm_matrix chi1inv_inf, eps_inf; + geps->eff_chi1inv_matrix(c, &chi1inv_inf, v, geps->tol, geps->maxeval, fallback, fill); + if (fallback) return false; + sym_matrix_invert(&eps_inf, &chi1inv_inf); + + std::complex tensor[9], tensor_inv[9]; + const double inf[9] = {eps_inf.m00, eps_inf.m01, eps_inf.m02, eps_inf.m01, eps_inf.m11, + eps_inf.m12, eps_inf.m02, eps_inf.m12, eps_inf.m22}; + for (int i = 0; i < 9; i++) + tensor[i] = std::complex(inf[i], 0); + + const struct { + material_type m; + double w; + } side[2] = {{front, fill}, {behind, 1.0 - fill}}; + + for (int sd = 0; sd < 2; sd++) { + if (!side[sd].m) continue; + for (const auto &su : side[sd].m->medium.E_susceptibilities) { + meep::lorentzian_susceptibility sus(su.frequency, su.gamma, su.drude); + // chi1 is linear in sigma, so evaluate the lineshape once at unit + // amplitude and scale. + std::complex lineshape = sus.chi1_discrete(freq, 1.0, geps->dt); + for (int i = 0; i < 9; i++) + tensor[i] += lineshape * side[sd].w * vec_to_value(su.sigma_diag, su.sigma_offdiag, i); + } + } + + vector3 zero = {0.0, 0.0, 0.0}; + for (int i = 0; i < 9; i++) { + double cond = 0; + for (int sd = 0; sd < 2; sd++) + if (side[sd].m) + cond += side[sd].w * vec_to_value(side[sd].m->medium.D_conductivity_diag, zero, i); + tensor[i] *= conductivity_factor(cond, freq, geps->dt); + } + + invert_tensor(tensor_inv, tensor); + int row = 0; + switch (component_direction(c)) { + case meep::X: + case meep::R: row = 0; break; + case meep::Y: + case meep::P: row = 1; break; + case meep::Z: row = 2; break; + case meep::NO_DIRECTION: return false; + } + for (int i = 0; i < 3; i++) + chi1inv_row[i] = tensor_inv[3 * row + i]; + return true; +} + void eff_chi1inv_row_disp(meep::component c, std::complex chi1inv_row[3], const meep::vec &r, double freq, geom_epsilon *geps) { std::complex tensor[9], tensor_inv[9]; @@ -3183,7 +3266,6 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf, std::vector fields_f, double *frequencies, double scalegrad, meep::grid_volume &gv, geom_epsilon *geps, int object_index, int *params, double du) { - (void)frequencies; (void)du; if (object_index < 0 || object_index >= geps->geometry.num_items) meep::abort("geometry_addgradient: object index %d out of range (%d objects)", object_index, @@ -3317,7 +3399,9 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf, double fill; const geometric_object *front = NULL; vector3 shiftby = {0, 0, 0}; - if (!geps->interface_fill(voxel, geps->tol, geps->maxeval, fill, &front, &shiftby)) + material_type mat_front = NULL, mat_behind = NULL; + if (!geps->interface_fill(voxel, geps->tol, geps->maxeval, fill, &front, &shiftby, + &mat_front, &mat_behind)) continue; if (front != obj) continue; if (fill <= 0.0 || fill >= 1.0) continue; @@ -3353,21 +3437,11 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf, for (int ax = 0; ax < 3; ax++) pixel_volume *= extent[ax]; - /* d(chi1inv)/d(fill), from meep's own tensor assembly. */ - bool fb_lo = false, fb_hi = false; - symm_matrix m_lo, m_hi; - geps->eff_chi1inv_matrix(adjoint_c, &m_lo, voxel, geps->tol, geps->maxeval, fb_lo, - std::max(0.0, fill - dfill)); - geps->eff_chi1inv_matrix(adjoint_c, &m_hi, voxel, geps->tol, geps->maxeval, fb_hi, - std::min(1.0, fill + dfill)); - if (fb_lo || fb_hi) continue; - - const double lo_row[3] = {m_lo.m00, m_lo.m01, m_lo.m02}; - const double hi_row[3] = {m_hi.m00, m_hi.m01, m_hi.m02}; - const double lo_row1[3] = {m_lo.m01, m_lo.m11, m_lo.m12}; - const double hi_row1[3] = {m_hi.m01, m_hi.m11, m_hi.m12}; - const double lo_row2[3] = {m_lo.m02, m_lo.m12, m_lo.m22}; - const double hi_row2[3] = {m_hi.m02, m_hi.m12, m_hi.m22}; + const double fill_lo = std::max(0.0, fill - dfill); + const double fill_hi = std::min(1.0, fill + dfill); + const double actual_dfill = fill_hi - fill_lo; + if (actual_dfill <= 0) continue; + int row_of = 0; switch (meep::component_direction(adjoint_c)) { case meep::X: @@ -3377,11 +3451,44 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf, case meep::Z: row_of = 2; break; default: continue; } - const double *lo_r = (row_of == 0) ? lo_row : (row_of == 1) ? lo_row1 : lo_row2; - const double *hi_r = (row_of == 0) ? hi_row : (row_of == 1) ? hi_row1 : hi_row2; - const double actual_dfill = std::min(1.0, fill + dfill) - std::max(0.0, fill - dfill); - if (actual_dfill <= 0) continue; - const double dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill; + + /* d(chi1inv)/d(fill), from meep's own tensor assembly. + + Which assembly depends on the materials at *this* interface, + not on md->trivial of the front object: a plain dielectric + sitting in front of a metal would otherwise take the + non-dispersive branch and drop the entire contribution at + exactly the pixels that carry the derivative. */ + std::complex dchi_dfill; + if (is_dispersive(mat_front) || is_dispersive(mat_behind)) { + std::complex lo_d[3], hi_d[3]; + if (!eff_chi1inv_row_disp_fill(adjoint_c, lo_d, voxel, frequencies[f_i], geps, + fill_lo, mat_front, mat_behind)) + continue; + if (!eff_chi1inv_row_disp_fill(adjoint_c, hi_d, voxel, frequencies[f_i], geps, + fill_hi, mat_front, mat_behind)) + continue; + dchi_dfill = (hi_d[dir_idx] - lo_d[dir_idx]) / actual_dfill; + } + else { + bool fb_lo = false, fb_hi = false; + symm_matrix m_lo, m_hi; + geps->eff_chi1inv_matrix(adjoint_c, &m_lo, voxel, geps->tol, geps->maxeval, fb_lo, + fill_lo); + geps->eff_chi1inv_matrix(adjoint_c, &m_hi, voxel, geps->tol, geps->maxeval, fb_hi, + fill_hi); + if (fb_lo || fb_hi) continue; + + const double lo_row[3] = {m_lo.m00, m_lo.m01, m_lo.m02}; + const double hi_row[3] = {m_hi.m00, m_hi.m01, m_hi.m02}; + const double lo_row1[3] = {m_lo.m01, m_lo.m11, m_lo.m12}; + const double hi_row1[3] = {m_hi.m01, m_hi.m11, m_hi.m12}; + const double lo_row2[3] = {m_lo.m02, m_lo.m12, m_lo.m22}; + const double hi_row2[3] = {m_hi.m02, m_hi.m12, m_hi.m22}; + const double *lo_r = (row_of == 0) ? lo_row : (row_of == 1) ? lo_row1 : lo_row2; + const double *hi_r = (row_of == 0) ? hi_row : (row_of == 1) ? hi_row1 : hi_row2; + dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill; + } const std::complex pair = std::complex(double(adj.real()), double(adj.imag())) * @@ -3400,7 +3507,7 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf, /* the leading minus matches get_material_gradient's convention, which returns -(d row/d parameter) */ local[nparams * f_i + ip] -= - node_weight * scalegrad * cyl_scale * dchi_dfill * dfill_dp * std::real(pair); + node_weight * scalegrad * cyl_scale * dfill_dp * std::real(dchi_dfill * pair); } } // node } diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index 5f485e9cd..ca1215f4c 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -233,7 +233,8 @@ class geom_epsilon : public meep::material_function { or fill < 0 if the pixel does not straddle an interface (so the shape derivative is zero there and the pixel can be skipped). */ bool interface_fill(const meep::volume &v, double tol, int maxeval, double &fill, - const geometric_object **which = NULL, vector3 *shift = NULL); + const geometric_object **which = NULL, vector3 *shift = NULL, + material_type *front_mat = NULL, material_type *behind_mat = NULL); void fallback_chi1inv_row(meep::component c, double chi1inv_row[3], const meep::volume &v, double tol, int maxeval);