Attempting 3D waveguide simulation

[Mythraptes]

Hello!

As part of my internship, I have been doing some work on MEEP, and we have gotten to the step where we are trying to gather the effective index of a waveguide with grating in 3D. Previously, we've done this with 2.5D (Sans grating) and received values that seemed to fit alternative method solutions fairly well. However, this is not the case for 3D.

Attempting to create a waveguide in 3D, then to use run_k_points to gather frequencies at certain points seems to strictly return the same range of frequencies for each k point, and not a consistent, incremental increase as is expected. Additionally, its especially difficult to get any frequencies below 0.5, and if I do, then I find a ceiling at 0.6. Through testing, I've found that, as far up the chain as I can understand it, it only reads reflected light? If I put a PML barrier in the z direction at the end of my cell, my k points return nothing despite being ahead of the barrier. Whereas if I dont put one, I get that same range of frequencies as mentioned previously.

I've brought in the code below. I'm practically at wit's end, I have tried nearly everything, tweaking nearly every parameter to get appropriate frequency values to come out, but whenever I do the math to obtain the effective index, the results are never right.

Does anyone have any idea how to get this to work?

import matplotlib.pyplot as plt
import numpy as np
import meep as mp
import math

freq_array = []
k_array = []
w_array = []

sy = 8  # size of cell in y direction
sx = 6
sz = 8

k_points = mp.interpolate(10, [mp.Vector3(0,0,0), mp.Vector3(0,0,0.5/sz)])
for k in k_points:
    k_array.append(k.z)

def main():
    # Some parameters to describe the geometry:
    eps = 13  # dielectric constant of waveguide
    w = 0.4  # width of waveguide
    n = 3.45
    n2 = 1.45
    lamb = 1.55
    alpha = 0.00004
    n_imag = lamb*alpha / (4*math.pi)

    #fcen = 1/lamb  # pulse center frequency
    fcen = 0.1
    df =  1.0 # pulse freq. width: large df = short impulse
    
    
    e_real = n**2 - n_imag**2 
    e_im = 2*n*n_imag
    d_cond = 2*math.pi * fcen * e_im / (e_real)

    e_real2 = n2**2 - n_imag**2
    e_im2 = 2*n2*n_imag
    d_cond2 = 2*math.pi*fcen*e_im2/e_real2
    

    # The cell dimensions

    dpml = 1  # PML thickness 

    cell = mp.Vector3(sx, sy, sz)

    b = mp.Block(size=mp.Vector3(0.4, 0.22, mp.inf),center=mp.Vector3(0, 0.11), material=mp.Medium(epsilon = e_real))#, D_conductivity = d_cond))
    d = mp.Block(size=mp.Vector3(mp.inf, 4, mp.inf),center=mp.Vector3(0, -2), material = mp.Medium(epsilon = n2**2))#, D_conductivity = d_cond2))
    


    
    s = mp.Source(
        src=mp.GaussianSource(fcen, fwidth=df),
        component=mp.Hz,
        center=mp.Vector3(0, 0, -3), #-1.885, and -1.95
    )

    #sym = mp.Mirror(direction=mp.X, phase=-1)

    

    sim = mp.Simulation(
        cell_size=cell,
        geometry=[b, d], #, c],
        sources=[s],
        #symmetries=[sym],
        boundary_layers=[mp.PML(dpml, direction = mp.X), mp.PML(dpml, direction = mp.Y)],#, mp.PML(dpml, direction = mp.Z)],
        resolution=10,
    )

    
    #sim.run(mp.after_sources(mp.Harminv(mp.Ez, mp.Vector3(), fcen, df)),
    #        until_after_sources=300)

    #sim.restart_fields()
    ##for k in k_points:
    print("Running simulation for k-points")
    freq_array.append(sim.run_k_points(100, k_points))  # Ensure t is specified
    


if __name__ == "__main__":
    main()

[stevengj]

Yes, 3d simulations work — we've done many, many validated 3d simulations. (But the more parameters you have, the easier it is to make a mistake in setting things up, and 3d simulations are harder to visualize and take longer to run so it's harder to debug them.)

I'm a little confused by your simulation. You are trying to get the dispersion relation of a waveguide? But then why wouldn't you use a 2d simulation, since your waveguide is invariant in one direction? i.e. why is sz > 0?

Also, it looks like you are exciting the wave with an Hz source at (0,0,-3), but then trying to detect the field by its Ez component at (0,0,0). This won't work because it looks like your structure has an x=0 mirror symmetry plane, in which case an Hz source is odd with respect to this plane and will excite only odd-symmetry fields, so your Ez component at x=0 will be zero … any nonzero Ez field you happen to detect will just be garbage from roundoff and perhaps discretization errors. In general, I tend to use harminv with the same field component and the same location as the source.

[Mythraptes]

Thank you for the response!

For why 3D with an invariant waveguide, this was intended as a test. We intend on graduating from the simple waveguide test to one with grating, and we can't do that with the pseudo-3d option.

And you're right that I made a mistake by trying to detect with Ez, just an artefact of having tried multiple things. I should have caught that, but it should be Hz for both. Also the symmetry was commented out, it was not in use at the time, I had just left it there on accident.

That said, even with the correction to Hz, the error remains the same. If we run this same simulation in pseudo-3d, then we get a set of frequencies that gradually increases with each k point (using run_k_points, thats the function that matters here) in a way we expect, and from there we are able to obtain sensible effective index values after doing the math. This doesn't happen with proper 3D, we get more or less the same set of values for each k point, and for the life of me I can't understand why.

[gummyequation]

I made a few changes to the code which helped me plot the frequencies and visualize the simulation. If you have any questions about the module functions I used, in the code let me know. You can get an idea of what the simulation looks like with sim.plot3D(). I used a few of my functions to plot the real parts of the frequencies and I am getting an incremental changes for the frequencies for every step-wise change in k-points. I do not know why you put the source outside the waveguide as stevengj has pointed out. But that said here's what I got when I widened to df = 3.0. As stevengj said, these can be nonsensical values I just want to help you visualize the simulation and manage the data better
(P.S. I left a lot of comments in your code as well as changes to help you move forward)

Here's what I got for df = 1, I may have increase nks,

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jul 19 15:42:14 2024

//Monday,July 22, 2024 Notes for githut_discussion_2865_copy.py
Testing for wider df values and lower res for quick results
@author: stinky
Where we are now jul 22, 15:11, want to test the current pulse with higher res

working to parallelize to speed up the process

"""
import matplotlib.pyplot as plt
import numpy as np
import meep as mp
import math
import hmfunc as hf # homemade functions for meep 
import pickle as pk
import os
import datetime
import time

# to get elapsed time for a sim
start = time.time()
timestamp = datetime.datetime.now().strftime("%Y-%m-%d-%-H:%M:%S")
print(timestamp)

freq_array = []
k_array = []
w_array = []

sy = 6  # size of cell in y direction
sx = 6
sz = 8


kband_pts = [mp.Vector3(0,0,0),mp.Vector3(0,0,0.5/sz)] # added this to for my plotting purposes 
nks = 25 # number of points between k-points
k_points= mp.interpolate(nks,kband_pts)# increased the number of points between the k-points


# for k in k_points:
#     k_array.append(k.z) I just use len(k_points)


# Some parameters to describe the geometry:
eps = 13  # dielectric constant of waveguide
w = 0.4  # width of waveguide
n = 3.45
n2 = 1.45
lamb = 1.55 
alpha = 0.00004
n_imag = lamb*alpha / (4*math.pi)


#fcen = 1/lamb  # pulse center frequency
fcen = 0.2
df = 3.0 # pulse freq. width: large df = short impulse


# Medium parameters
e_real = n**2 - n_imag**2 
e_im = 2*n*n_imag
d_cond = 2*math.pi * fcen * e_im / (e_real)

e_real2 = n2**2 - n_imag**2
e_im2 = 2*n2*n_imag
d_cond2 = 2*math.pi*fcen*e_im2/e_real2



# The cell dimensions
dpml = 1  # PML thickness 
res = 10 # resolution
cell = mp.Vector3(sx, sy, sz)

b = mp.Block(size=mp.Vector3(0.4, 0.22, mp.inf),center=mp.Vector3(0, 0.11), material=mp.Medium(epsilon = e_real))#, D_conductivity = d_cond))
d = mp.Block(size=mp.Vector3(mp.inf, 4, mp.inf),center=mp.Vector3(0, -2), material = mp.Medium(epsilon = n2**2))#, D_conductivity = d_cond2))

# Soruces
s = mp.Source(
    src=mp.GaussianSource(fcen, fwidth=df),component=mp.Ez, 
    center=mp.Vector3(0, 0, -1),) #-1.885, and -1.95

#sym = mp.Mirror(direction=mp.X, phase=-1)


sim = mp.Simulation(
    cell_size=cell,
    geometry=[b, d], #, c],
    sources=[s],
    k_point =k_points, # specificies to the simulation that there are periodic boundary conditions
    #symmetries=[sym],
    boundary_layers=[mp.PML(dpml, direction = mp.X), mp.PML(dpml, direction = mp.Y)],#, mp.PML(dpml, direction = mp.Z)],
    resolution=res,)


# sim.plot3D()  # 3-D render of simulation

#sim.run(mp.after_sources(mp.Harminv(mp.Ez, mp.Vector3(), fcen, df)), # TE
#        until_after_sources=300)
# sim.restart_fields()
#for k in k_points:
    
# print("Running simulation for k-points")
all_freqs= sim.run_k_points(100, k_points) # instead of loop I just make a single reference, for a single array 
real_freqs = hf.get_real_parts(all_freqs) #just splits the array from the imag parts and keeps the real


# pickle export
# useful variables to export when analyzing the pickled data in a seperate script, this way you do not have to run the simulation everytime
export_variables = {
    'a' : sz,
    "kband_pts": kband_pts,
    "nks" : nks,
    "cell" : cell,
    "resolution" : res,
    "fcen": fcen,
    "df": df,
    "dpml": dpml,
    "k_points": k_points,
    }

file_path= os.path.abspath(__file__) #absolute path

# Extract the filename without the path
filename = os.path.basename(file_path)
filename_no_type = filename.removesuffix('.py')+timestamp+"_res_" +str(res)
juicy_pickle = f'{filename_no_type}.pickle'

# outputs the all_freqs data into binary 
with open(juicy_pickle, "wb") as f: # wb = write binary 
    pk.dump((all_freqs,export_variables), f) # f is the placeholder from above with open(...) as "f"
    
# plotting 
# hf.plot_freqs(real_freqs,kband_pts) # external function 

end=time.time()

time_elapsed_sec_min = '\nelapsed run time: '+ f'{end-start}' +' seconds'+ "\nElaspsed time in minutes: "+ f"{(end-start)/60}"
print(time_elapsed_sec_min)

   time_elapsed_sec_min = '\nelapsed run time: '+ f'{end-start}' +' seconds'+ "\nElaspsed time in minutes: "+ f"{(end-start)/60}"
    print(time_elapsed_sec_min)

[Mythraptes]

Thank you for the response!

I greatly appreciate the detail you put into this, unfortunately the problem persists. Speaking of the source, it actually does not matter whether or not it is located in or out of the waveguide, the values will somehow still be the same, at each k point. That is part of what is slowly driving me insane.

I've attached a graph of what we expect the frequency vs k points graph to look like (we've obtained this through pseudo-3D simulations) and we know this data is close to correct as we've validated it elsewhere. If we adapt that same pseudo-3D code to actual 3D, then the problem reappears.