Yee Lattice of meep at the boundary of two media

[Baobabs-kuang]

Hello there friends from the meep community,
We have been working on a PINN(Physics Informed Neural Network) recently, and we need to substitute the electric field obtained from the solution in meep into Maxwell equations to calculate the LOSS function, but I found that when I substitute the EM field data and dielectric constant data derived from meep calculations into the equations, it doesn't conform very well 1-2 pixels on the boundary of the two media, as shown in the following two diagrams, the first one is the structure diagram of my test structure diagram, where the green box is the PML, the red line is the source, and the blue box is the region where we verified Maxwell equations, and the second figure is the residuals of importing the EMF data derived in meep into the curl equation, which is close to 0 in most regions, but the residuals are larger at the boundary of the media
I get the EMF and dielectric constant information by calling sim.get_dft_array and sim.get_array
This question has been bothering me for many days and I hope I can get a hand from you.

[oskooi]

The features in your "residual" image are likely due to artifacts of the bilinear interpolation Meep uses internally to center the DFT fields of the array slice on the Yee-grid voxels.

(#2625 updates the docstrings for the array-slice routines to make this property explicit.)

[Baobabs-kuang]

Thanks for the reply, but I think you may have misunderstood my question.
What I mean is that meep solves Maxwell's equations using FDTD, so the DFT field and dielectric constant derived by meep should match Maxwell's equations, but when I use get_array and get_dft_array to derive the dielectric constant and DFT field data to substitute into the curl equation, it matches better in the non-boundary dielectric region, but in the two or more dielectric boundary 1-2 pixel region, the exported data do not conform to the curl equation! My residual definition is as shown in the figure below, where the superscript with meep refers to the dielectric constant and DFT field data derived through get_array and get_dft_array in meep.

[Baobabs-kuang]

I guess the dielectric constant data derived through get_array is not the data used by meep when calculating the difference equation? Otherwise, in any case, the DFT field and dielectric constant data calculated by meep are substituted into the curl equation, and the equation should be established!
I'm not quite sure if my guess is reasonable, and I'm wondering how to derive the correct dielectric constant data

[Baobabs-kuang]

Here's the code for simulation.

import numpy as np
import meep as mp
from scipy.sparse import diags
import matplotlib.pyplot as plt

def Grident(array, dim):
    if dim == 0:
        grident_array_lens =array.shape[0]
        grident_array = diags([1, -1, 1, -1], [-grident_array_lens+1, -1, 1, grident_array_lens-1],
                              shape=(grident_array_lens, grident_array_lens))/2
        gradient = grident_array@array
        return -gradient
    if dim == 1:
        grident_array_lens = array.shape[1]
        grident_array = diags([1, -1, 1, -1], [-grident_array_lens + 1, -1, 1, grident_array_lens - 1],
                              shape=(grident_array_lens, grident_array_lens))/2
        gradient = np.transpose(grident_array @ np.transpose(array))
        return gradient

def cal_DFT(eps_array, res):
    gp = 0.1
    gh = 0.5
    dpml = 1
    dsub = 0.5
    dpad = 0.5
    lcen = 0.98
    fcen = 1/lcen
    df = 0.2*fcen

    SiO2 = mp.Medium(index=1.5)

    sx = gp*len(eps_array)
    sy = dpml+dsub+gh+dpad+dpml

    resolution = res

    k_point = mp.Vector3(0, 0, 0)
    pml_layers = [mp.PML(thickness=dpml, direction=mp.Y)]
    cell_size = mp.Vector3(sx, sy, 0)

    src_pt = mp.Vector3(0, -0.5*sy+dpml+0.5*dsub, 0)
    mon_pt = mp.Vector3(0, 0.5 * sy - dpml - 0.5 * dpad, 0)

    geometry = [mp.Block(material=SiO2, size=mp.Vector3(mp.inf, dpml+dsub, mp.inf), center=mp.Vector3(0, -0.5*sy+0.5*(dpml+dsub), 0))]
    source = [mp.Source(mp.GaussianSource(fcen, fwidth=df, is_integrated=True), component=mp.Ez, center=src_pt, size=mp.Vector3(sx, 0, 0))]

    for i in range(len(eps_array)):
        geometry.append(
            mp.Block(material=mp.Medium(index=eps_array[i]),
                     size=mp.Vector3(gp, gh, 0),
                     center=mp.Vector3(-sx/2+gp/2+i*gp, -0.5*sy+dpml+dsub+gh/2))
        )

    sim = mp.Simulation(
        resolution=resolution,
        cell_size=cell_size,
        boundary_layers=pml_layers,
        geometry=geometry,
        k_point=k_point,
        sources=source
    )

    calculate_arrange = mp.Volume(center=mp.Vector3(), size=mp.Vector3(sx, sy - 2 * dpml - dsub*1.2))
    dft_obj = sim.add_dft_fields([mp.Hx, mp.Hy, mp.Ez], fcen, 0, 1, where=calculate_arrange)

    sim.run(until_after_sources=mp.stop_when_fields_decayed(3, mp.Ez, mon_pt, 1e-5))

    Ez = sim.get_dft_array(dft_obj, mp.Ez, 0)
    Hx = sim.get_dft_array(dft_obj, mp.Hx, 0)
    Hy = sim.get_dft_array(dft_obj, mp.Hy, 0)
    epsilon = sim.get_array(component=mp.Dielectric, vol=calculate_arrange)

    field = np.empty((6, *Ez.shape))
    field[0, ...] = Ez.real
    field[1, ...] = Ez.imag
    field[2, ...] = Hx.real
    field[3, ...] = Hx.imag
    field[4, ...] = Hy.real
    field[5, ...] = Hy.imag

    return field, epsilon


eps_array = np.ones(10)
eps_array[0:2] = 2

res = 50
field, epsilon = cal_DFT(eps_array, res)
omega = 2*np.pi/0.98

DxHy_real = res*Grident(field[4, 1:-1, :], 0)
DxHy_imag = res*1j*Grident(field[5, 1:-1, :], 0)
DyHx_real = res*(Grident(field[2, 1:-1, :], 1))
DyHx_imag = res*1j*Grident(field[3, 1:-1, :], 1)
Ez_real = field[0, 1:-1, :]
Ez_imag = 1j*field[1, 1:-1, :]

residual = DxHy_real[:, 1:-1] + DxHy_imag[:, 1:-1] + DyHx_real[:, 1:-1] + DyHx_imag[:, 1:-1] - 1j * omega * epsilon[1:-1, 1:-1] * (Ez_real[:, 1:-1]+Ez_imag[:, 1:-1])
residual_abs = np.abs(residual)


plt.figure()
plt.pcolor(residual_abs)
plt.colorbar()
plt.clim(0,5)
plt.title('residual')
plt.show()

[stevengj]

I guess the dielectric constant data derived through get_array is not the data used by meep when calculating the difference equation?

Right. get_array with mp.Dielectric returns a kind of average scalar dielectric constant at the center of each voxel, which is much simpler than the actual dielectric constant used internally. In particular, Meep does anisotropic subpixel averaging of ε near material interfaces, as described in this paper, so internally its ε is a tensor, with different tensor components (from different averages) at every Yee lattice point.

Moreover, when you get fields with sim.get_dft_array(dft_obj, mp.Ez, …), it is also interpolated from the Yee lattice to the center of each voxel. Though there is an option yee_grid=True to compute the fields on the underlying Yee grid(s) without interpolation.

So, the upshot is that neither the fields nor the ε you are getting are quite what Meep uses internally for the discretized curl equations, due to the interpolation. (Basically, the reason for all of this interpolation is that raw Yee grids are quite tricky for users to deal with.)

[Baobabs-kuang]

thanks!!