Meep

目的

安装完Meep之后一直很少用它来做具体的计算,因为COMSOL和Lumerical FDTD Solutions确实有图形界面非常好用,最近想要计算偶极子耦合到直波导的耦合效率,用COMSOL计算三维结构需要较大的内存,用FDTD的时候服务器没有图形界面,扫描设置很不方便,这时候决定好好打磨一下Meep以为我所用。我的思路是:

  1. 学会用Meep计算光波导模式并和COMSOL对上
  2. 学会用Meep设置偶极子的源计算并计算器耦合到波导模式的功率

今天笔记主要讲的是第一步。

那里可以学习?

Meep如果直接去浏览官方文档还是比较低效的,最好是可以直接从例子中学习,有不懂再去看。我也希望我把我的例子写出来可以帮助大家快速上手,用Meep计算波导模式有两个方法,我是各种浏览github的问题区域知道的,里面别人分享了很多的源代码。比如别人问为啥MPB和Meep计算的结果不一样,我就可以知道有两个方法计算波导模式:

总之就是多看别人的例子再去查询函数用法,我觉得是最快的。

实现过程

我现在要计算的是,宽度为650nm,高度为250nm的SiN波导在波长为780nm处的波导模式。其结构如下

MPB实现

首先是插入一下包和设置基本的结构参数

1
2
3
4
5
6
7
8
9
10
import meep as mp
from meep import mpb
import numpy as np
from matplotlib import pyplot as plt

# Then the wavelength
omega = 1/0.78
material_SiN=mp.Medium(index=2.02)
width_guide=0.65
height_guide=0.25

接下来就是要初始化处几何结构(geometry),格点个数(resolution,每1um多少个),我们先用MPB实现,这样我们还需要初始化三维原胞的结构,如下

1
2
3
4
5
6
7
8
resolution=200
size_cal_y=2
size_cal_z=2

geometry = [mp.Block(material=material_SiN,
size=mp.Vector3(mp.inf,width_guide, height_guide),
center=mp.Vector3(0, 0, height_guide/2))]
geometry_lattice = mp.Lattice(size=mp.Vector3(0, size_cal_y, size_cal_z))

我们的波导方向沿着X,只是为了计算模式,厚度设置为0,为一个二维问题。然后是设置我们要计算几个模式以及初始化计算对象
1
2
3
4
5
6
7
8
num_modes = 4

ms = mpb.ModeSolver(
geometry_lattice = geometry_lattice,
geometry = geometry,
resolution = resolution,
num_bands = num_modes
)

参数的详细设置可以参见源文档。接下来我们需要设置一个函数来讲每次计算的模式的场分布取出,这是通过上面的链接issues851学到的,即需要设置
1
2
E = []
store_fields = lambda ms_temp, mode: E.append(ms_temp.get_efield(which_band=mode,bloch_phase=False))

然后就可以计算了
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
k = ms.find_k(
mp.NO_PARITY,
omega, # omega
1, # band_min
num_modes, # band_max
mp.Vector3(1,0,0), # korig_and_kdir
1e-4, # tol
omega * 2, # kmag_guess
omega * 0.1, # kmag_min
omega * 2, # kmag_max
store_fields, # band_funcs
)
neff=k[0]/omega # 有效折射率
eps = ms.get_epsilon() # 折射率分布
eps_arr=np.transpose(np.array(eps))

绘制结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# %%
x_list=np.linspace(-size_cal_y/2,size_cal_y/2,resolution*2)
y_list=np.linspace(-size_cal_z/2,size_cal_z/2,resolution*2)
x_grid,y_grid=np.meshgrid(x_list,y_list)

plt.figure(figsize=(12,4))
# Plot the E fields
for mode in range(num_modes):
Ex=np.squeeze(E[mode][:,:,0,1]).transpose()
print('Current band: {}'.format(mode+1))
plt.subplot(1,num_modes,1+mode)
plt.pcolormesh(x_grid,y_grid,eps.transpose(), cmap='binary')
plt.pcolormesh(x_grid,y_grid,np.abs(Ex), cmap='jet', alpha=0.9)
plt.axis('off')
st_title='Ex of mode '+str(mode+1)+' from MPB neff=' +str('%.2f' %(k[mode]/omega))
plt.title(st_title)
plt.tight_layout()
plt.savefig('SimulatedFieldMPB.png')
plt.show()


上面分别绘制出来了四个模式的Ex分量的实部。

Meep实现

首先是插入一下包和设置基本的结构参数,和之前一样

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import meep as mp
from meep import mpb
import numpy as np
from matplotlib import pyplot as plt

#Then the wavelength
omega = 1/0.78
material_SiN=mp.Medium(index=2.02)
width_guide=0.65
height_guide=0.25
resolution=200
size_cal_y=2
size_cal_z=2

geometry = [mp.Block(material=material_SiN,
size=mp.Vector3(mp.inf,width_guide, height_guide),
center=mp.Vector3(0, 0, height_guide/2))]

不同的是接下来我们需要设置PML边界条件,并同样的初始化对象
1
2
3
4
5
6
7
8
9
10
11
12
13
num_modes = 4

dpml = 0.5
pml_layers = [mp.PML(dpml)]
sim = mp.Simulation(cell_size=mp.Vector3(0,size_cal_y+dpml*2,size_cal_z+dpml*2),
resolution=resolution,
geometry=geometry,
boundary_layers=pml_layers,)

sim.init_sim()
where = mp.Volume(center=mp.Vector3(), size=mp.Vector3(0,size_cal_y,size_cal_y))
direction = mp.X #波导方向
kpoint = mp.Vector3(1,0,0)

然后是计算,但是得提前准备好坐标,因为meep的函数get_eigenmode函数只能得到一个amplitude的方法,你可以通过该方法插值出场的数值,这个也是我仔细研究了一下函数的用法才发现的,所以算模式还是MPB吧
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
x_list=np.linspace(-size_cal_y/2,size_cal_y/2,resolution*2)
y_list=np.linspace(-size_cal_z/2,size_cal_z/2,resolution*2)
x_grid,y_grid=np.meshgrid(x_list,y_list)
pos_extract=mp.Vector3(np.ravel(x_grid)*0,np.ravel(x_grid),np.ravel(y_grid))
Ex_mp=np.zeros((resolution*2,resolution*2,num_modes),dtype=complex)
# Extract the field in a loop
num_modes=4
k_meep=np.zeros(num_modes)
for mode in range(num_modes):
print('Simultion Mode'+str(mode+1))
ed = sim.get_eigenmode(omega, direction, where, mode+1, kpoint, eig_vol=None, match_frequency=True,
parity=mp.NO_PARITY, resolution=resolution, eigensolver_tol=1e-12)
k_meep[mode]=ed.k[0]
for l in range(resolution*2):
for m in range(resolution*2):
Ex_mp[l,m,mode]=ed.amplitude(mp.Vector3(0,x_list[m],y_list[l]),1)

然后绘制结果
1
2
3
4
5
6
7
8
9
10
11
plt.figure(figsize=(14,4))
for mode in range(num_modes):
print('Current band: {}'.format(mode+1))
plt.subplot(1,num_modes,1+mode)
plt.pcolormesh(x_grid,y_grid,np.abs(Ex_mp[:,:,mode]), cmap='jet', alpha=0.9)
plt.axis('off')
st_title='Ex of mode '+str(mode+1)+' from MEEP neff=' +str('%.2f' %(k_meep[mode]/omega))
plt.title(st_title)
plt.tight_layout()
plt.savefig('SimulatedFieldMeep.png')
plt.show()

放在一起比较

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
plt.figure(figsize=(12,8))
for mode in range(num_modes):
print('Current band: {}'.format(mode+1))
plt.subplot(2,num_modes,1+mode)
plt.pcolormesh(x_grid,y_grid,np.abs(Ex_mp[:,:,mode]), cmap='jet', alpha=0.9)
plt.axis('off')
st_title='Ex of mode '+str(mode+1)+' from MEEP neff=' +str('%.2f' %(k_meep[mode]/omega))
plt.title(st_title)
for mode in range(num_modes):
Ex=np.squeeze(E[mode][:,:,0,1]).transpose()
print('Current band: {}'.format(mode+1))
plt.subplot(2,num_modes,1+mode+num_modes)
plt.pcolormesh(x_grid,y_grid,np.abs(Ex), cmap='jet', alpha=0.9)
plt.axis('off')
st_title='Ex of mode '+str(mode+1)+' from MPB neff=' +str('%.2f' %(k[mode]/omega))
plt.title(st_title)

plt.tight_layout()
plt.savefig('SimulatedFieldCompare.png')
plt.show()

可以看见二者结果一样,我也与COMSOL比对,结果是对的上的。因此通过Meep计算波导模式还是很简单的,这是一个最简答的波导模式,下一个笔记将会介绍如何计算更加复杂结果的波导模式以及需要注意的问题,以及计算偶极子耦合效率。

源代码

可以直接看我的github
WaveguideModeSimulations

觉得有帮助的话记得点赞加关注吧~