A dispersive MaterialGrid is silently non-dispersive in the forward simulation

[smartalecH]

A MaterialGrid whose constituents carry E_susceptibilities behaves as if it had none. The poles are interpolated by epsilon_material_grid and then never reach the timestepping.

Reproduced on 1.35.0-beta. A block of a medium against a MaterialGrid(air, medium) with all weights at 1, which should be indistinguishable:

	Drude	index 1.5 (control)
block of the medium	3.76353163e-02	6.86560106e+00
MaterialGrid at u=1	1.66692970e+00	6.86560106e+00
vacuum	1.66692970e+00	1.66692970e+00

The non-dispersive control matches to all printed digits, so the projection and the instantaneous-epsilon interpolation are fine. The dispersive case matches vacuum, bit for bit.

import numpy as np, meep as mp

DRUDE = mp.Medium(epsilon=1.0, E_susceptibilities=[
    mp.DrudeSusceptibility(frequency=1.0, gamma=0.5, sigma=1.0)])

def transmission(material):
    sim = mp.Simulation(
        cell_size=mp.Vector3(6, 4), resolution=20,
        boundary_layers=[mp.PML(1.0)],
        sources=[mp.Source(mp.GaussianSource(1.0, fwidth=0.2), component=mp.Ez,
                           center=mp.Vector3(-1.8, 0))],
        geometry=[mp.Block(center=mp.Vector3(), size=mp.Vector3(1, 1),
                           material=material)],
        eps_averaging=False, force_complex_fields=True)
    mon = sim.add_dft_fields([mp.Ez], 1.0, 0, 1,
                             center=mp.Vector3(1.8, 0), size=mp.Vector3())
    sim.run(until=150)
    return float(np.abs(sim.get_dft_array(mon, mp.Ez, 0)) ** 2)

grid = mp.MaterialGrid(mp.Vector3(2, 2), mp.air, DRUDE, do_averaging=False, beta=0)
grid.update_weights(np.ones(4))

print(transmission(DRUDE), transmission(grid), transmission(mp.air))

Cause

geom_epsilon::add_susceptibilities registers poles only for materials that satisfy is_medium:

for (int i = 0; i < geometry.num_items; ++i)
  if (is_medium(geometry.items[i].material, &mm))
    pols = add_pols(pols, ft == meep::E_stuff ? mm->E_susceptibilities : ...);

and is_medium returns true only for which_subclass == MEDIUM. A material grid is MATERIAL_GRID, so nothing ever calls sigma_row for its poles and the amplitudes epsilon_material_grid interpolates go nowhere. No other path registers them; extra_materials is documented for material_function, not for materials attached to geometric objects.

D_conductivity is unaffected, because geom_epsilon::conductivity reads md->medium through get_material_pt rather than going through pole registration. That is probably why this has gone unnoticed.

Consequence for the adjoint solver

get_chi1_tensor_disp reads medium.E_susceptibilities directly through get_material_pt, so it does see the poles. Optimizing a dispersive material grid therefore returns a confident gradient for a structure that was never simulated: a central difference of the objective with respect to a design variable comes back exactly zero while the adjoint reports a large value.

Related

While reproducing this, mu, H_susceptibilities, the nonlinearities and B_conductivity also turn out never to be written to the mixed medium — epsilon_material_grid interpolates only epsilon, the electric susceptibilities and the diagonal D_conductivity. Those are accepted and then ignored, which is the same failure mode.

Fix in #3286.