reproducing slow-light waveguide in MPB

[stevengj]

Hi @fengwen1206, I'm currently trying to reproduce your slow-light waveguide results in MPB, and it's not quite matching.

In particular, this is ε(x,y) in my (periodic) computational cell in MPB:

loaded from your design file and mapped to a refractive index of n_Si = 3.476.

which gives me a band diagram (Hz polarization) of

This looks qualitatively similar to what you have, but quantitatively there is a mismatch of 2–3% in the dispersion relation of the guided band

which is bigger than I naively expected for a resolution of 40 pixels/a:

Not sure what could explain the difference.

On the other hand, for the Hz polarizations the discretization errors arising from dielectric boundaries tend to be larger, due to the discontinuity in the E field for this polarization, but I can't really go any higher since you only provided data at 40 pixels/a.

Can we try a cross-check at a higher resolution?

[fengwen1206]

Hi Steven,

I actually use the square elements, will this cause the big frequency shift?

[stevengj]

I'm using a completely different discretization, not FEM at all — MPB is based on a planewave basis. Of course, different schemes should converge to one another as the resolution increases. It's arguably a better check to use completely different discretizations than to try another code (e.g. Comsol) using an identical discretization, since the latter does not help to quantify the discretization error.

(A basic difficulty is that the discontinuity in the E field tends to reduce the rate of convergence of every scheme to 1st-order unless interfacial discontinuities are specifically taken into account).

[stevengj]

If I plot all of the bands, you'll see that the discrepancies get larger for higher bands (= smaller λ), which is consistent with the problem being discretization error:

[fengwen1206]

How can we solve the problem? Cross-check with high resolutions? If so, how would you like me to generate a higher-resolution design?

I could extract the contour plot of the design and simulate it using body-fitted mesh. Will it work better for MPB?

[stevengj]

Yes, if you extract the contour plot of the design, then you can use that to generate an ε(x,y) array (e.g. in a CSV file) at a much higher resolution (e.g. 4x resolution) which I can then use in MPB (and hopefully get 1/4 the error since I will be seeing 1st-order convergence), and you can simulate with your own code (at a higher resolution and/or a fitted mesh). This will also give us an estimate of the discretization error in the original optimization.

(PS. The list-of-(x,y,value) format is somewhat inconvenient to deal with in MPB or Meep. Ideally I would just like a 2d array of values on a rectangular grid, as a 2d CSV file. No xy coordinates are needed since that can be inferred from the computational cell size.)

[stevengj]

(I tried doubling the resolution in MPB with your existing design file — just linearly interpolating — and it doesn't change my results significantly. So it's really the precise description of the boundary location that will matter.)

[stevengj]

I computed the difference more quantitatively, and it's actually only 2–3% for the slow-light band, which isn't too bad:

[fengwen1206]

I have updated the design as array format and also upload the design pattern in a high resolution. I have created pull request for the updates.

As I used square elements, height/width=408/40=10.2 which is a bit less than 6 sqrt(3). I imported the contour of the design into comsol and compared with comsol for the guided mode for ka/2pi=0.4, the relative error is around 1.5%. Comsol's result is also slightly smaller than my result (40 pixels/a). I will get the full band information tomorrow morning.

Slightly different volume fraction of Si in the supercell due to discretization will cause frequency shift. Dilated design will lead to a slight lower frequency and eroded design will results in higher frequency.

[stevengj]

When I run MPB at a resolution of 80 pixels/a using your high-resolution image, with a cell width of 10.2 (instead of 6√3), I get a relative error of around 1.5% (the MPB frequencies are systematically lower than yours):

[stevengj]

I ran MPB at a resolution of 160 pixels/a and the numbers changed by only about 0.004%, so it seems pretty well converged at 80 pixels/a.

[stevengj]

In contrast, if I increase the cell with to 6√3 instead of 10.2 (stretching the whole structure by 1.9%), then the difference with your code increases to about 2.5%, very similar to my previous results using your low-resolution design.

So, the difference in cell width explains nearly half of the discrepancy, and the remaining 1.5% difference is attributable to discretization error.

[fengwen1206]

I did resolution convergence study compared with Comsol with a constant cell width of 10.2. I used the high resolution as input and reduce the resolution from 160 pixels/a to 80 pixels/a and 40 pixels/a by averaging (Fengwen-ReEval, blue curves). The interpolation scheme is 1/ε_e=1/ε_air+x_e ( 1/ ε_si -1/ε_air) as in the optimization.

[fengwen1206]

following the same approach, if I use the interpolation scheme is ε_e= ε_air+x_e ( ε_si - ε_air) ( Fengwen-ReEval, blue curves). I obtained the result below with very small errors compared to comsol result. The design pattern in160 pixels/a is black-white and the result is same as above.

In the optimization procedure, I use interpolation in 1/ε. So when I extract the contour line, I actual need to extract contour line in ε. I will update the designs and corresponding results later.

I have attached the band data from Comsol. Your result should match very well with it.
Comsol_Band.csv

[fengwen1206]

I tried to re-extract contour design so that the averaged permittivity is very close between the optimized one and the extract contour design. Redo all the comsol and matlab evaluatoion. Using the interpolation scheme is ε_e= ε_air+x_e ( ε_si - ε_air) ( Fengwen-ReEval, blue curves), I get the results very closed to the optimized result both for re-evaluation in my code and comsol. As expected, the relative error for group index is slightly higher than frequency.

The new high resolution design is
HigRes_DesMatch_Opt_Dnum_2.csv

Band strucuture from comsol evaluation is
Comsol_Band_HigRes_DesMatch.csv
Group index from comsol evaluation is
Comsol_Group_HigRes_DesMatch.csv

[fengwen1206]

Shall we also include the high resolution design in the manuscript?

[stevengj]

I think we should just mention something about the cross-check/validation in the text, say something about the expected accuracy at a given resolution, and include the high-resolution design an an MPB validation script in the github repository?

[fengwen1206]

Sound good. I think we should use the second high resolution design, Could you please evaluate this design?
Thanks.
HigRes_DesMatch_Opt_Dnum_2.csv

[stevengj]

Re-ran using HigRes_DesMatch_Opt_Dnum_2.csv. The difference is now about 0.3% (using sx = 10.2).

---

My MPB code, for reference:

import numpy as np
import meep as mp
from meep import mpb
import matplotlib.pyplot as plt
import scipy.interpolate as interp
import h5py

sx = 10.2 # not quite 6*np.sqrt(3)
# sx = 6*np.sqrt(3)

eps_hi = np.loadtxt("HigRes_DesMatch_Opt_Dnum_2.csv", delimiter=',')
eps_hi.shape

plt.imshow(eps_hi)
plt.colorbar()
plt.axis('off')
plt.show()

hf = h5py.File('HigRes_DesMatch_Opt_Dnum_2.h5', 'w')
hf.create_dataset('eps', data=eps_hi.T)
hf.close()

ms = mpb.ModeSolver(num_bands=18,
                    eigensolver_block_size=12,
                    k_points=mp.interpolate(20, [mp.Vector3(0,0), mp.Vector3(0,0.5)]),
                    epsilon_input_file='HigRes_DesMatch_Opt_Dnum_2.h5',
                    geometry_lattice=mp.Lattice(size=mp.Vector3(sx, 1)),
                    resolution=80)

ms.init_params(mp.TE, True)
eps2 = ms.get_epsilon()
plt.imshow(eps2.T, cmap='binary')
plt.axis('off')
plt.colorbar()
plt.show()

ms.run_te()

design_results = np.loadtxt("Opt_Band.csv", delimiter=',')

ky = [k.y for k in ms.k_points]
te_freqs = ms.all_freqs
plt.plot(ky, te_freqs, "g-");
plt.title("band structure (MPB = green, Fengwen = red)")
plt.plot(design_results[:,0], design_results[:,1:], "r-");
plt.xlabel("$ka/2\\pi$")
plt.ylabel("$\\omega a/2\\pi c$")
plt.legend(["Fengwen", "MPB"])
plt.show()

plt.plot(ky, te_freqs[:,12], "go-");
plt.plot(design_results[:,0], design_results[:,13], "r*--");
plt.title("guided bands: MPB vs Fengwen")
plt.xlabel("$ka/2\\pi$")
plt.ylabel("$\\omega a/2\\pi c$")
plt.legend(["Fengwen", "MPB"])
plt.show()

hz = np.real(ms.get_hfield(13)[:,:,0,2])
plt.imshow(hz.T, cmap="RdBu")
plt.show()

design_interp = interp.interp1d(design_results[:,0], design_results[:,1:], axis=0)(ky)
plt.title("% error: MPB (resolution 80 pixel/a) vs Fengwen")
plt.ylabel("$|\\Delta\\omega|/\\omega$ (%)")
plt.xlabel("$ka/2\\pi$")
plt.plot(ky, 100*(np.abs(design_interp - te_freqs) / te_freqs)[:,12], '.-');
plt.show()