Simulation results don't match up with HFSS simulation results

[10UAikawa]

Hi MEEP community, I have tried to resolve my problem of getting the correct transmission results for my periodic structure. I have used my structure that I have in HFSS to verify my results in MEEP. However, the results I am getting in MEEP do not seem to be giving me the correct results like they do in HFSS.

I am not quite sure where my error is arising from, but I have a few ideas.

-I have defined my materials in the THz range as instructed in this part of the documentation: https://meep.readthedocs.io/en/latest/Materials/#conductivity-and-complex

-When I remove the materials conductive materials the results remain the same as if they were not there. (It should be noted I am working with thin films so the thickness of the ribbon in comparison to the rest of the structure could maybe be an issue, idk)

-Maybe the way in which I am doing normalization might be incorrect. I normalize my structures results in the second simulation with respect to the Gaussian transmitting through a vacuum. (I had an earlier MEEP Q&A tell me this is the correct way though)

-Maybe the way in which I am defining the periodicity might not be correct. My structure should only be periodic in the x and y direction.

-Lastly doing the simulation takes a much short time to do in HFSS compared to MEEP. I am currently considering 1um as my MEEP unit. Could I change my MEEP unit to cm and thus be able to increase the resolution more? Because the size of my unit cell ends up being 400um x 400um x 2400um which leads to something that is very computationally expensive.

Here is what my structure looks like:

Here are the results that I should get for my transmission:

Here are the results I am getting with and without the light blue ribbon (my conductive material in the THz range) in MEEP:

Any suggestions on my possible problem will help immensely!

Thanks in Advance!

Here is the Code:

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

'''-----units are in microns (1 MEEP unit = 1 um)-----'''
'''-----with time scale as t=3.3e-15s since a/c gives time scale-----'''
wavelength = 749.481 # wavelength in microns for a frequency at 0.4THz

'''-------varibles to change that effect simulation time---------'''
resolution = 0.1 # pixels per mircon (pixel/um)
decay_by = 1e-3 # decay from the maximum feild
nfreq = 500 # number of frequencies points at which to compute flux
dT = 100000 # keeps running this additional time units once cirteria has been meet
'''--------------------------------------------------------------'''

#Padding
dpml = 50

#size of the cell
X_size = 400
Y_size = 400
Z_size = 800+dpml*2+wavelength*2

#dimensions of the quartz block
x_block = 400
y_block = 400
z_block = 100

#dimensions of the ribbon 
ribbon_width = 230
rib_thickness = 0.02

#dimensions of the aluminum
alum_slit_spacing = 344
alum_slit_width = 5   
alum_slit_height = 220
alum_thickness = 0.02

#Materials
Cd3As2_conductivity = 1e5  # S/m
Cd3As2_epsilon = 1  # relative permitvity is already dimensionless (1 for metals/semimetals)
Al_conductivity = 38000000 # S/m
Al_epsilon = 1.000021 # relative permitvity is already dimensionless (1 for metals/semimetals)
Quartz_epsilon = 4.41  # relative permitvity is already dimensionless 

#frequnecy definitions
wave_max = 428.275 # wavelength of frequency at 0.7THz
wave_min = 1498.962 # wavelength of frequency at 0.2THz
freq_intrest = 1/wavelength # The frequency of the light source is calculated as the reciprocal of the wavelength, given that the speed of light in Meep units is one.
f_max = 1/wave_max # wavelength of frequency at 0.7THz [MEEP units]
f_min = 1/wave_min # wavelength of frequency at 0.2THz [MEEP units]
df = f_max-f_min #range of frequencies

# conversion = σD=(a/c)σ/(εrε0) [from MEEP docmentation]
eo=8.8541878128E-12  #in F/m [SI units]
a=1e-6  #units in meters from micron to meters (since we use micron units as 1) [SI units]
c=(3e8) #in (m/s) [SI units]

# Using a constant conductivty because we are examining this a narrow frequency band so we can use this as an approximtion
Cd3As2_conductivity_MEEP_units = (a/c*Cd3As2_conductivity/(Cd3As2_epsilon*eo))
Al_conductivity_MEEP_units = (a/c)*Al_conductivity/(Al_epsilon*eo)

# imaginary part of the permitivty
Eps_img_Cd3As2 = Cd3As2_epsilon*Cd3As2_conductivity_MEEP_units/(freq_intrest*(2*np.pi))
Eps_img_Al = Al_epsilon*Al_conductivity_MEEP_units/(freq_intrest*(2*np.pi))

# Create a material with the specified conductivity/permitivty
Cd3As2_material = mp.Medium(epsilon = Cd3As2_epsilon, D_conductivity = (2*np.pi)*freq_intrest*Eps_img_Cd3As2/Cd3As2_epsilon) 
Al_material = mp.Medium(epsilon = Al_epsilon, D_conductivity = (2*np.pi)*freq_intrest*Eps_img_Al/Al_epsilon)  
Quartz_material = mp.Medium(epsilon = Quartz_epsilon)
Vacuum_material = mp.Medium(epsilon=1)

#defining cell size
cell = mp.Vector3(X_size, Y_size, Z_size)

#Bloch-periodic boundary conditions
k_point = mp.Vector3(x_block/2, y_block/2, 0) #because periodic in the x and y direction 

#boundary conditions
pml_layers = [mp.Absorber(thickness=dpml,direction=mp.Z)] # because periodic boundary conditions

# locations for the transmittance and reflectance                   ORDER of the positions
source_loc = -(Z_size/2)+dpml*5 # source location                   TOP
reflection_loc = -(Z_size/2)+dpml*6 # reflectance location          FIRST MIDDLE
transmittance_loc = (Z_size/2)+dpml*1 # transmittance location      SECOND MIDDLE
pt_loc = transmittance_loc = (Z_size/2)-dpml*2 # pt location        BOTTOM

#Normailzing WRT the incident plane wave
geometry = [
    #setting to a vaccum
    mp.Block(
            mp.Vector3(X_size, Y_size, Z_size),
            center=mp.Vector3(0,0,0),
            material=Vacuum_material,
        )

]

#setting up the source 
sources = [
    mp.Source(
        mp.GaussianSource(freq_intrest, fwidth=df),
        component=mp.Ez,
        center=mp.Vector3(0,0, source_loc ),  
        size=mp.Vector3(X_size, Y_size, 0),
    )
]

#setting up the simulation
sim = mp.Simulation(
    cell_size=cell,
    boundary_layers=pml_layers,
    geometry=geometry,
    k_point=k_point,
    sources=sources,
    resolution=resolution,
)

# reflected flux
refl_fr = mp.FluxRegion(center=mp.Vector3(0,0,reflection_loc), size=mp.Vector3(x_block,y_block,0))

# specifying the location of the monitor for the reflectance
refl = sim.add_flux(freq_intrest, df, nfreq, refl_fr)

# transmitted flux
tran_fr = mp.FluxRegion(center=mp.Vector3(0,0,transmittance_loc), size=mp.Vector3(x_block,y_block,0))

# specifying the location of the monitor for the transmittance
tran = sim.add_flux(freq_intrest, df, nfreq, tran_fr)

# the point at which the feilds must decay sufficiently
pt = mp.Vector3(0,0,pt_loc) 

# running simulation until the feild decays sufficiently
sim.run(until_after_sources=mp.stop_when_fields_decayed(dT,mp.Ez,pt,decay_by=decay_by))

# for normalization run, save flux fields data for reflection plane
straight_refl_data = sim.get_flux_data(refl)

# save incident power for transmission plane
straight_tran_flux = mp.get_fluxes(tran)

'''---reset to include structure---'''
sim.reset_meep()

geometry = [
    #Quartz (bottom)
    mp.Block(
        mp.Vector3(x_block, y_block, z_block),
        center=mp.Vector3(0,0,0),
        material=Quartz_material,
    ),

    # # Cd3As2 Ribbon (top)
    # mp.Block(
    #     mp.Vector3(ribbon_width, y_block, rib_thickness),
    #     center=mp.Vector3(0,0,(z_block/2)+(rib_thickness/2)),
    #     material=Cd3As2_material,
    # ),   
]

#setting up the source 
sources = [
    mp.Source(
        mp.GaussianSource(freq_intrest, fwidth=df),
        component=mp.Ez,
        center=mp.Vector3(0,0, source_loc ),  
        size=mp.Vector3(X_size, Y_size, 0),
    )
]

#setting up the simulation
sim = mp.Simulation(
    cell_size=cell,
    boundary_layers=pml_layers,
    geometry=geometry,
    k_point=k_point,
    sources=sources,
    resolution=resolution,
)

# reflected flux
refl = sim.add_flux(freq_intrest, df, nfreq, refl_fr)

# transmitted flux  (no need to redefine since same location)
tran = sim.add_flux(freq_intrest, df, nfreq, tran_fr)

# for normal run, load negated fields to subtract incident from refl. fields
sim.load_minus_flux_data(refl, straight_refl_data)

# the point at which the feilds must decay sufficiently
pt = mp.Vector3(0,0,pt_loc)

# running simulation until the feild decays sufficiently
sim.run(until_after_sources=mp.stop_when_fields_decayed(dT,mp.Ez,pt,decay_by=decay_by))

# gather fluxs from the monitors
refl_flux = mp.get_fluxes(refl)
tran_flux = mp.get_fluxes(tran)
flux_freqs = mp.get_flux_freqs(refl)

wl = []
Rs = []
Ts = []

# Filling arrays with data and normailzing results WRT to the vaccum (i.e when Epsilon=1)
for i in range(nfreq):
    wl = np.append(wl, (3e3)/(1/flux_freqs[i]) ) # converts to THz
    Rs = np.append(Rs,-refl_flux[i]/straight_tran_flux[i])
    Ts = np.append(Ts,tran_flux[i]/straight_tran_flux[i])

# Plotting reflectance, transmittances, and loss
if mp.am_master():
    plt.figure()
    plt.plot(wl,Rs,'bo-',label='reflectance')
    plt.plot(wl,Ts,'ro-',label='transmittance')
    plt.plot(wl,1-Rs-Ts,'go-',label='loss')
    plt.xlabel("Frequency (THz)")
    plt.legend(loc="upper right")
    plt.show()

# Getting e-feild data 
ez_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Ey)
# Choose the z-slice (e.g., z_index) you want to visualize
index_slice = 0
# Extract the 2D slice in the x-y plane at the chosen z-slice
ez_slice = np.real(ez_data[:, index_slice,:])


'''-------------This is how I can view a slice of a 3d object in a 2d plane------------------'''
# Define the slice plane
slice_pos = (z_block/2)+(rib_thickness)+(z_block)+(alum_thickness)  # Position of the slice along the z-axis
slice_dim = mp.Vector3(0, Y_size, Z_size) # y-z plane

# Get the Dielectric data (for example, the Dielectric component to see the geometry)
Die_data = sim.get_array(center=mp.Vector3(0, 0, slice_pos), size=slice_dim, component=mp.Dielectric)

# Plot the 2D slice using Matplotlib
plt.figure()
plt.imshow(Die_data.transpose(), interpolation='spline36', cmap='binary')
plt.imshow(ez_slice.transpose(), interpolation='spline36', cmap='RdBu', alpha=0.7)
plt.axis('on')
plt.show()
'''------------------------------------------------------------------------------------------'''

[manujmg]

I am dealing with the same problem. Your problem resides in the resolution you are using.
Your resolution is 0.1 pixels/um which is around 150 wavelengths (quite large). Therefore, the minimum dimension you can use for any block would be 1/resolution which in your case would be 1/0.1(pixels/um) = 10um/pixel. Anything smaller than that will not be considered in your simulation.

Try to change the thicknesses of those materials to 1/resolution so they are considered in your simulation.

This is exactly what I am doing right now to represent PEC materials with a thickness lower than 1/resolution. My question here for MEEP developers would be....

1. Is it possible to assign different resolutions for each coordinate?

• resolution_x_y = 10 #for x,y axis
• resolution_z = 1000 #for z axis

2. Is it possible to assign different resolutions for different structures? Let's say I have a geometry made of three blocks, can I change the mesh (resolution) of the third block?

I guess one workaround is to download the source code and modify the mesh generation code. But, that will require a decent amount of time.

[10UAikawa]

Thanks for the reply. Do you by chance get this warning when initializing your structure: "Warning: grid volume is not an integer number of pixels; cell size will be rounded to nearest pixel."

[peterropac]

I am not a developer, but I know the answers to your questions. People have asked them before...

1. No, having different resolutions in different directions is impossible.
2. no, it is impossible to have different resolutions in different regions, at least not in the same simulation.

[manujmg]

Yes, see below
"
Initializing structure...
Warning: grid volume is not an integer number of pixels; cell size will be rounded to nearest pixel.
Padding y to even number of grid points.
Halving computational cell along direction y
time for choose_chunkdivision = 0.000364065 s
Working in 3D dimensions.
"