WIP: mode-solver improvements #333
Conversation
…ks for both sources and mode decomposition
…memory leak in mode coefficients
|
Is this ready to be merged? We have a test involving a binary grating (below) for which the analytic result is known using scalar diffraction theory. The diffraction/transmission efficiency is [2/(πq)]^2 for q odd where q is the diffraction order (the q even modes have zero efficiency). Thus, as an example, the row with "grating:, " from the output with "1" (for q) in the second column should have a fourth column of [2/(π*1)] ^2 = 0.40528, etc. import meep as mp
import math
resolution = 20 # pixels/um
nglass = 1.5
glass = mp.Medium(index=nglass)
dsub = 3 # substrate thickness
dair = 3 # air padding between grating and pml
gp = 10.0 # grating period
gh = 0.5 # grating height
gdc = 0.5 # grating duty cycle
lambda_min = 0.4 # vacuum wavelength
lambda_max = 0.8 # vacuum wavelength
fmin = 1/lambda_max
fmax = 1/lambda_min
fcen = 0.5*(fmin+fmax)
df = fmax-fmin
dpml = lambda_max # PML thickness
sx = dpml+dsub+gh+dair+dpml
sy = gp
cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
geometry = [ mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)),
mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0)) ]
k_point = mp.Vector3(0,0,0)
src_pos = -0.5*sx+dpml
sources = [ mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=mp.Vector3(src_pos,0,0), size=mp.Vector3(0,sy,0)) ]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources)
# fmode is the frequency corresponding to free-space wavelength lambda=0.535 um
fmode = 1.0/0.535
xm = 0.5*sx-dpml
mflux = sim.add_eigenmode(fmode, 0, 1, mp.FluxRegion(center=mp.Vector3(xm,0,0), size=mp.Vector3(0,sy,0)))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(xm,0), 1e-6))
##################################################
# kpoint_function to compute propagation vector
# for diffraction order #mode at frequency freq
##################################################
def kpoint_func(freq, mode):
k0 = freq
ky = (mode-1)/gp
if ky<k0:
kx=math.sqrt(k0*k0-ky*ky)
return mp.Vector3(kx,ky,0.0)
return mp.Vector3(0.0,0.0,0.0)
##################################################
# get eigenmode coefficients of propagating modes
##################################################
for nm in range(1,10):
kpoint = kpoint_func(fmode, nm)
if kpoint.x > 0.0:
alpha = sim.get_eigenmode_coefficients(mflux, [nm], kpoint_func)
alpha0Plus = alpha[2*0 + 0] # coefficient of forward-traveling fundamental mode
alpha0Minus = alpha[2*0 + 1] # coefficient of backward-traveling fundamental mode
angle=math.atan2(kpoint.y,kpoint.x) * 180/math.pi
if (mp.am_master()):
print("grating:, {}, {}, {}, {}, {}".format(fmode,nm-1,angle,abs(alpha0Plus**2),abs(alpha0Minus**2))) |
|
The following is the output of the above test (3rd, 5th, and 6th columns only which show the diffraction order and mode coefficients for the forward and backward propagating fundamental mode) for the current master branch and this PR. Current Master Branch (latest commit: "silence a few integer conversion warnings") PR #333 There is a slight difference but nevertheless the updated results still do not agree with the scalar diffraction theory. An even simpler check is that the last column should all be zeros. |
|
You should be using #294 is also a problem here, since you want to really only look for TM modes — there are twice as many modes as you think, because the eigenmode coefficients are computing both polarizations here. Also, there is another factor of two in the number of modes because you are forgetting that each diffraction order is doubled (one going in the +y direction and one in the -y direction). For both these reasons, your Also, I'm not sure what scalar-diffraction method you are using, but I have a hard time believing this is an exact result — your supposedly exact answer doesn't even depend upon the refractive index of the glass. (If you set |
|
I've pushed a fix for #294. You can now pass an |
|
(Note, by the way, that the documentation for |
|
The example is from Prof. George Barbastathis' 2.71 Optics class on OCW: see slide 9 of the lecture 16 notes. I was able to reproduce these results using the near2far approach as shown in these notes. The right figure on the first slide shows the relative amplitude of the diffraction peaks which agrees with the scalar theory. |
|
@ChristopherHogan, once this PR is merged we'll probably want to update mode_coeffs.py to a more stringent test — set |
|
@oskooi, did you have the wrong link? The slides you linked don't have anything about diffracted orders, and slide 9 is about ray tracing. You're not going to convince me that there is an exact answer independent of |
|
Oops, wrong link. Corrected now. |
|
By the way, I don't see why |
|
To use the result from the 2.71 notes, you have to choose the grating height S so that the phase difference is π between propagating a distance S in the glass vs. a distance S in the air. i.e. the thickness S must depend on both the index of the glass and the wavelength. It doesn't look like your code does this. (Even if you choose S correctly to get the π binary phase, this is only an approximation, because it neglects diffraction within the thickness S of the grating, which is only accurate in the limit where the pitch Λ is much larger than both S and the wavelength λ. It doesn't look like you are operating in this limit, so I don't think the approximation will be super accurate for you even with the correct S.) |
|
The pitch Λ is 10 μm, the wavelength λ is 0.535 μm, and the grating height S is 0.5 μm so we are well within the scalar limit. |
|
I reran the test above using the latest version of this PR (commit "whoops") with slight modification to the last few lines to add the for nm in range(1,10):
alpha = sim.get_eigenmode_coefficients(mflux, [nm], eig_parity=mp.TM)
alpha0Plus = alpha[2*0 + 0] # coefficient of forward-traveling fundamental mode
alpha0Minus = alpha[2*0 + 1] # coefficient of backward-traveling fundamental mode
if (mp.am_master()):
print("grating:, {}, {}, {}, {}".format(fmode,nm-1,abs(alpha0Plus**2),abs(alpha0Minus**2)))The output, showing just the last three columns, is: The scalar theory is valid for λ=0.5 μm (which produces a π phase difference in the grating of height S=0.5 and index 1.5 with air). The results for the λ=0.535 μm, which is what we are currently using, should be similar. Obviously, the results are still incorrect. Rather than focus on the second column which is the square of the mode coefficient for the forward-propagating fundamental mode (i.e., the diffraction efficiency), it would be better to find out why the last column (the squared mode coefficient of the backward-propagating fundamental model) is not 0. One thing to notice in the MPB's Why, as you claimed earlier, should ky be 0? Also, note that the number of iterations for convergence is not a small number. |
|
ky=0 is equivalent to m/gp: these two differ by a reciprocal lattice vector. In any case, ky is ignored because it is set by the boundary conditions. (The fact that all valid ky values are equivalent to 0 is the basis of diffraction theory: the only allowed diffracted orders are the planewave modes that are compatible with the periodic boundary conditions. Equivalently, in a periodic system, the ky of the incident planewave is conserved modulo reciprocal lattice vectors.) |
|
The coefficient of the backward-propagating mode won't be zero because of finite resolution. If you increase the resolution it should decrease. (This is for the same reason that an eigenmode source does not produce purely a forward-propagating mode at a finite resolution.) |
|
I ran the test at several resolutions (pixels/μm): 20, 40, 80, 160, 200. The backward mode coefficients do not seem to be changing noticeably. The forward mode coefficients are not converging to the expected result. resolution=20 resolution=200 How can we verify that MPB is computing the correct mode since since we cannot use the wavevector from |
Then maybe they are limited by something other than resolution, e.g. the PML thickness or the finite runtime or the eigensolver tolerance (which I think defaults to 1e-7). The coefficients are small enough that they look like they must be some sort of approximation error.
How can you tell? To compare to the (2/πq)² answer you need to:
Why can't you look at kₓ? |
|
The wavevectors from MPB's The first five kx values should therefore be: The values reported by MPB (at a resolution of 200 in Meep) are quite different: Thus, it seems that there is some mixup in the modes MPB is calculating with the ones it should be calculating for this analysis. Perhaps this is related to adding Another indication that MPB is computing the wrong modes is that the mode coefficients of the forward-propagating mode for the even diffraction orders should all be zero. The output from earlier shows that these mode coefficients are not even close to zero. |
|
@oskooi, as I said above, for each kx there are actually two modes, at ±ky, which is why the MPB solutions are in pairs. When you account for this, the MPB modes are matching the analytic results to 1e-4. Hence, it appears MPB is computing exactly the correct modes. |
|
Two modes at ±ky for each kx would suggest that the "band indices" parameter passed to The mode coefficients of the forward-propagating mode are quite different for the +ky and -ky diffraction orders even though the planewave is normally incident; how come? Given the symmetry of the problem, shouldn't they be equal? Also, how can we reduce MPB's eigensolver tolerance and increase resolution to improve the accuracy of the results? Currently, there is no way to set these parameters through |
It's exactly the band number in MPB, which refers to the modes in order of ω. As you say, this is not the diffraction order m, and it can't be because the diffraction order is signed (corresponding to ±ky).
Because of the degeneracy, MPB is giving you a "random" superposition of the ±ky solutions. The sum of their powers should be the same as the sum of the +ky and -ky powers. If you pass the EVEN_Y parity, it will only compute the +ky + -ky (cosine) mode. (For ODD_Y, it will compute the sine mode, but this will have zero power because your source is even.) At ky≠0 the degeneracy breaks and then each mode is a planewave. You could also put in a ky that is slightly ≠ 0, say 1e-6, to get planewave modes again. (Remember, you change ky via the boundary conditions.) |
Yes, it should really have all of the same parameters (including eig_resolution and eigensolver_tol) as add_eigenmode_source. |
|
Yikes, I just noticed that |
|
I created a new def kpoint_func(freq, nmode):
k0 = freq
if nmode == 1:
return mp.Vector3(k0,0,0)
else:
[d, m] = divmod(nmode,2)
ky = d/gp
kx=math.sqrt(k0*k0-ky*ky)
return mp.Vector3(kx,ky+1e-6,0)
for nm in range(1,8):
alpha = sim.get_eigenmode_coefficients(mflux, [nm], eig_parity=mp.ODD_Z, kpoint_func=kpoint_func)
alpha0Plus = alpha[2*0 + 0] # coefficient of forward-traveling fundamental mode
alpha0Minus = alpha[2*0 + 1] # coefficient of backward-traveling fundamental mode
if (mp.am_master()):
if nm == 1:
print("grating:, {}, {}, {}".format(nm,abs(alpha0Plus**2),abs(alpha0Minus**2)))
else:
[d, pm] = divmod(nm,2)
print("grating:, {}, {}, {}, {}, {}".format(nm,d,pm,abs(alpha0Plus**2),abs(alpha0Minus**2)))Even when taking into account the transmissivity of the air/glass interface and multiplying by cos(θ) = kₓc/ω for each diffracted mode, the results (at resolution=80) are still not in agreement: Modes 4 (+ky) and 5 (-ky) for diffraction order m=2 should both have zero coefficient for the forward-propagating mode. Mode 4 seems to be close but not 5. Also, the output from MPB is still showing ky=0 although this may not be an indication that the degeneracy persists: It would be good to rerun this test by reducing MPB's tolerance and/or increasing resolution. Would it be possible to add a commit to this PR to make these parameters as keyword arguments (and also set their default values to be consistent with those in |
* mpb match-freq should only search k in direction d * adjust for bloch-periodic boundaries in get_eigenmode, so that it works for both sources and mode decomposition * allow get_eigenmode to return NULL if mode was not found, clean up a memory leak in mode coefficients * whoops * add parity argument to get_eigenmode_coefficients (fises NanoComp#294) * add missing get_eigenmode_coefficients parameters, fix defaults
| if (match_frequency && kmatch == 0) { | ||
| vec cen = eig_vol.center(); | ||
| kmatch = omega_src * sqrt(get_eps(cen)*get_mu(cen)); | ||
| k[d-X] = kmatch * R[d-X][d-X]; // convert to reciprocal basis |
There was a problem hiding this comment.
This line seems problematic and may be the cause of the incorrect mode reported in #2107.
The fix: kmatch should be changed to (kmatch - kpara_sum) where kpara_sum is the sum of the magnitude of the parallel components of k (i.e., parallel to the monitor). Instead, this line assigns a single component (in the normal direction d) of k to be kmatch (in Cartesian coordinates). The resulting k therefore has the wrong magnitude if any of its other components are non-zero. This may explain why the mode found by MPB is not what is expected.
There was a problem hiding this comment.
It's just an initial guess for the Newton iteration; it doesn't have to be that accurate.
No description provided.