Simulation results change with resolution

[mawc2019]

It seems that simulation results may not quickly converge as resolution increases. Here is an example, where the transmission spectrum is calculated as the normalized transmission power. The central wavelength is 0.5. The structure is as follows.

The transmission spectra for a series of resolution are as follows. The unit of resolution is the number of points per unit length.

The normalized transmission power at the central frequency (=2) under a series of resolution is as follows. If we assume that the result at the resolution 200 is correct, we need the resolution to be higher than 90 in order to attain a relative error below 10%, and we need the resolution to be higher than 110 in order to attain a relative error below 5%. Is this situation normal?

The code for calculating the transmission spectrum at a given resolution is as follows.

import meep as mp
import numpy as np
from matplotlib import pyplot as plt
mp.verbosity(0)
Si,SiO2 = mp.Medium(index=3.4),mp.Medium(index=1.44)

resolution = 40
Sx,Sy = 1,8
block_width,monitor_height = Sx,2
pml_layers = [mp.PML(1.0,direction=mp.Y)]
cell_size = mp.Vector3(Sx,Sy)

fcen,fwidth,nf = 2,0.4,200
source_center,source_size  = [0,-2,0],mp.Vector3(Sx,0)
src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.EigenModeSource(src,eig_band = 1,size = source_size,center=source_center)]

geometry = [mp.Block(center=mp.Vector3(), size=mp.Vector3(block_width,block_width), material=Si)]
sim = mp.Simulation(cell_size=cell_size,boundary_layers=pml_layers,geometry=geometry,sources=source,k_point=mp.Vector3(),
                    default_material=SiO2,resolution=resolution)

fluxregion = mp.FluxRegion(center=mp.Vector3(0,monitor_height),size=mp.Vector3(Sx,0))
trans = sim.add_flux(fcen,fwidth,nf,fluxregion)

fig = plt.figure()
sim.plot2D(ax=fig.gca())
plt.show()

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(0,monitor_height), 1e-8))
power = np.array(mp.get_fluxes(trans))
sim.reset_meep()

sim = mp.Simulation(cell_size=cell_size,boundary_layers=pml_layers,geometry=[],sources=source,k_point=mp.Vector3(),
                    default_material=SiO2,resolution=resolution)
fluxregion = mp.FluxRegion(center=mp.Vector3(0,monitor_height),size=mp.Vector3(Sx,0))
trans = sim.add_flux(fcen,fwidth,nf,fluxregion)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(0,monitor_height), 1e-8))
power_ref = np.array(mp.get_fluxes(trans))
sim.reset_meep()

freqs = np.linspace(fcen-fwidth/2,fcen+fwidth/2,nf)
plt.plot(freqs,power/power_ref)
plt.title("Transmission spectra")
plt.xlabel("Frequency")
plt.ylabel("Normalized transmission power");

[MarkMa1990]

according to numerical stability of FDTD, you have to have like wavelength/20/index, here, e.g., 500nm/20/3.4, almost resolution close to 135. this should always be in mind.

[mawc2019]

Thanks! But the issue remains.

you have to have like wavelength/20/index, here, e.g., 500nm/20/3.4, almost resolution close to 135. this should always be in mind.

Indeed, satisfying such a stringent requirement is safe, but it does not reveal the rate of convergence as the resolution increases. I wonder whether the rate of convergence is normal. To be specific, in the test, when the resolution amounts to 10 points per wavelength/index=0.5/3.4, the relative error is above 25%. Is such a large error at such a level of resolution normal? By the way, in some built-in examples in Meep, they do not use the resolution as high as 20 points per wavelength/index. For example, the resolution here amounts to 9.2 points per wavelength/index=1.56/3.4.

according to numerical stability of FDTD

This issue does not seem a matter of numerical stability. In FDTD, numerical stability typically requires that the time step has an upper bound related to the space increment in the lattice. In Meep, as the default value of the Courant number is 0.5, the time step is always adjusted accordingly unless it is otherwise specified.

[MarkMa1990]

Thanks! But the issue remains.

you have to have like wavelength/20/index, here, e.g., 500nm/20/3.4, almost resolution close to 135. this should always be in mind.

Indeed, satisfying such a stringent requirement is safe, but it does not reveal the rate of convergence as the resolution increases. I wonder whether the rate of convergence is normal. To be specific, in the test, when the resolution amounts to 10 points per wavelength/index=0.5/3.4, the relative error is above 25%. Is such a large error at such a level of resolution normal? By the way, in some built-in examples in Meep, they do not use the resolution as high as 20 points per wavelength/index. For example, the resolution here amounts to 9.2 points per wavelength/index=1.56/3.4.

according to numerical stability of FDTD

This issue does not seem a matter of numerical stability. In FDTD, numerical stability typically requires that the time step has an upper bound related to the space increment in the lattice. In Meep, as the default value of the Courant number is 0.5, the time step is always adjusted accordingly unless it is otherwise specified.

Oh, sorry, It is called the numerical dispersion that requires space discretization as dx (in 3D) < lambda / 12 / index ~ close to resolution 82. in experience, the lambda here should be the minimum wavelength of the light source (at max frequency in gaussian source). so, in your case, lambda is around 454nm, the required resolution should be at least 90. in fact, about the time sampling, dt < T/12 should be fulfilled as well. So, to me, FDTD has two conditions to make good simulations: one is the courant conditon dt< 1/sqrt(3) dx; the other is space discretization, dx < lambda (in medium) / 12. as soon as all these two conditions are fulfilled, you can get close to the true propagations. otherwise, the light phase speed will mess up and the simulation can not be trusted.

[mawc2019]

It is called the numerical dispersion that requires space discretization as dx (in 3D) < lambda / 12 / index ~ close to resolution 82. in experience, the lambda here should be the minimum wavelength of the light source (at max frequency in gaussian source). so, in your case, lambda is around 454nm, the required resolution should be at least 90.

Although the source contains a range of frequencies, the two black curves in the last figure above (where I just supplemented the curve for the relative error) were evaluated at a single frequency, which was the central frequency. When the range of frequency is shrunk from 0.4 to 0.1, the change in these two curves are negligible.

one is the courant conditon dt< 1/sqrt(3) dx; the other is space discretization, dx < lambda (in medium) / 12. as soon as all these two conditions are fulfilled, you can get close to the true propagations.

Sometimes we may not obtain satisfactory results when the two conditions are fulfilled. In the previous example, the relative error is about 15% when the two conditions are fulfilled (and the time step obeys a more stringent condition by default). It does not seem acceptable.
In some cases the situation may be even worse. In the following example, the sum of adjoint gradients with respect to the design region is calculated, which should converge as resolution increases. The simulation involves an EigenmodeCoefficient adjoint solver and was run on v1.23.0 (because recent versions may have some bugs that produces incorrect adjoint gradients). The geometric structure and media are the same as before. The results are shown below. As shown by the blue curves for the sum of adjoint gradients, even when the resolution is 140, which amounts to 20.6 points per wavelength/index, the relative error can be as high as -30%.

The code for calculating the sum of adjoint gradients using an EigenmodeCoefficient adjoint solver at a given resolution is as follows. This example does not involve extra normalization. When the normalization is included in the same way as the previous example, the relative errors do not change noticeably.

import meep as mp
import meep.adjoint as mpa
import autograd.numpy as npa
import numpy as np

Si,SiO2 = mp.Medium(index=3.4),mp.Medium(index=1.44)

resolution,design_region_resolution = 40,40

Sx,Sy = 1,8
design_region_width,design_region_height = 1,1
pml_layers = [mp.PML(1.0,direction=mp.Y)]

cell_size = mp.Vector3(Sx,Sy)

fcen,fwidth = 1/0.5,0.4
source_center,source_size  = [0,-2,0],mp.Vector3(Sx,0)
src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.EigenModeSource(src,eig_band = 1,size = source_size,center=source_center)]

Nx,Ny = int(design_region_resolution*design_region_width),int(design_region_resolution*design_region_height)

design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')
design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_region_width, design_region_height, 0)))

geometry = [mp.Block(center=design_region.center, size=design_region.size, material=design_variables)]

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

emc = mpa.EigenmodeCoefficient(sim,mp.Volume(center=mp.Vector3(0,2,0),size=mp.Vector3(x=Sx)),mode=1)
ob_list = [emc]

def J(emc):
    power = npa.abs(emc) ** 2
    return power

opt = mpa.OptimizationProblem(simulation = sim,objective_functions = J,objective_arguments = ob_list,
                              design_regions = [design_region],fcen=fcen,df = 0,nf = 1,decay_by=1e-10)

x0 = 1.0*np.ones(Nx*Ny)
opt.update_design([x0])

f0, dJ_dn = opt()
sim.reset_meep()

dJ_sum = np.sum(dJ_dn)
print("The sum of adjoint gradients is ",dJ_sum)

[MarkMa1990]

when the two conditions are fulfilled, in theory, it guarantees error within around 1~2%. However, this numerical error accumulates while propagating, so resolutions might need to be higher to get that accuracy. In terms of gradients, I think there are some problems are not solved, especially for bloch-boundary where the manually introduced boundary condition break the Lorentz-transformation. However, while doing reversing k-point in simulation settings or k-point are all zero, it should work. based on what I know of meep, it has problems as follows:
#2054
#2087

I have a feeling, it is the setting of adjoint source that introduce the problem, which seems to be a big issue. I tested for 3D a few months ago, and it has big error as well. I need to find out that result.

[stevengj]

I'm not sure there is actually an issue here — the key point is that when you have a phenomenon involving a sharp resonance, even a tiny shift in the resonant frequency can cause a huge change in the transmission. In your plot at the top, at a resolution of 40 (about 20 pixels/λ in air but only about 7 pixels/λ in the dielectric), the resonant frequency is shifted by about 10%, which is actually pretty good, but of course this causes a huge change in the transmitted power at any particular frequency.

This also happens in practice, because even if your Meep calculation is infinitely precise, when you go to manufacture a structure the resonant frequencies will inevitably be shifted by a few percent. To compensate, the experiment needs to either be able to tune their operating frequency or to tune the structure (e.g. thermally).

[stevengj]

So, in practice, what you are interested in is the convergence of the peak value.

Above, for a resolution of 40 (= 7 pixels/λ in the dielectric), your resonant frequency is off by about 10% and the peak transmission is also off by about 10%. This seems quite reasonable to me.

[stevengj]

As for adjoint gradients, you don't want the adjoint gradient to match the infinite resolution (exact Maxwell) solve. You want the gradient of your (approximate) forward solve, so that optimization can evolve the parameters to improve the forward solve.