Replies: 2 comments 1 reply
-
Here is my codeimport meep as mp
import meep.adjoint as mpa
import autograd.numpy as npa
from autograd import tensor_jacobian_product, grad
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.patches import Circle
import nlopt
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
seed = 240
np.random.seed(seed)
mp.quiet(quietval=True)
InP = mp.Medium(index=3.1882)
InP_PQ = mp.Medium(index=3.3555)
Air = mp.air
##As before, we'll define our geometry, filtering, and wavelength parameters.
waveguide_width = 2
design_region_width = 20
design_region_height = 4
arm_separation = 1.0
waveguide_length = 3
pml_size = 3.0
resolution = 20
t_InP = 4 #bottom
t_InP_PQ = 0.4 #waveguide
design_region_resolution = int(resolution)
frequencies = 1 / np.linspace(1.31, 1.61, 5)
##We'll also define our simulation domain and set up a broadband source.
Sx = 2 * pml_size + 2 * waveguide_length + design_region_width
Sy = Sx/2
Sz = 2 * pml_size + t_InP_PQ + 4
cell_size = mp.Vector3(Sx, Sy, Sz)
pml_layers = [mp.PML(pml_size)]
fcen = 0.5*(1 / 1.31 + 1 / 1.61)
fwidth = 1.31 - 1 / 1.61
source_center = [-Sx / 2 + pml_size + waveguide_length / 3 +2, 0, 0]
source_size = mp.Vector3(0, 3, 3)
kpoint = mp.Vector3(1, 0, 0)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
mp.EigenModeSource(
src,
eig_band=1,
direction=mp.NO_DIRECTION,
eig_kpoint=kpoint,
size=source_size,
center=source_center,
)
]
##Next, we'll define our design region. This time, however, we'll include a symmetry across the Y plane (normal direction of the symmetry plane points in the Y direction).
Nx = int(design_region_resolution * design_region_width) + 1
Ny = int(design_region_resolution * design_region_height) + 1
design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny), Air, InP_PQ, grid_type="U_MEAN", do_averaging = False) # changed
design_region = mpa.DesignRegion(
design_variables,
volume=mp.Volume(
center=mp.Vector3(),
size=mp.Vector3(design_region_width, design_region_height, t_InP_PQ), #changed
),
)
##We also need to define a bitmask that describes the boundary conditions of the waveguide and cladding layers.
##Define spatial arrays used to generate bit masks
x_g = np.linspace(-design_region_width / 2, design_region_width / 2, Nx)
y_g = np.linspace(-design_region_height / 2, design_region_height / 2, Ny)
X_g, Y_g = np.meshgrid(x_g, y_g, sparse=True, indexing="ij")
##Define the core mask
left_wg_mask = (X_g == -design_region_width / 2) & (np.abs(Y_g) <= waveguide_width / 2) #INPUT WAVEHUIDE
top_wg_mask = (Y_g == design_region_width / 2) & (np.abs(X_g) <= waveguide_width / 2) #OUTPUT WAVEHUIDE
Si_mask = left_wg_mask | top_wg_mask
##Define the cladding mask
border_mask = (
(X_g == -design_region_width / 2)
| (X_g == design_region_width / 2)
| (Y_g == -design_region_height / 2)
| (Y_g == design_region_height / 2)
)
SiO2_mask = border_mask.copy()
SiO2_mask[Si_mask] = False
##Finally we can formulate our geometry and simulation object.
input_wg = mp.Block(
center = mp.Vector3(x=-Sx / 4),
material = InP_PQ,
size = mp.Vector3(Sx / 2 + 1, waveguide_width, t_InP_PQ)
)
output_wg = mp.Block(
center = mp.Vector3(x= Sx / 4),
material = InP_PQ,
size = mp.Vector3(Sx / 2 + 1, waveguide_width, t_InP_PQ)
)
design_R = mp.Block(
center=design_region.center, size=design_region.size, material=design_variables
)
Bottom = mp.Block(
center=mp.Vector3(z = -t_InP_PQ/2-t_InP/2),
material=InP,
size=mp.Vector3(Sx, Sy, t_InP)
)
geometry = [
input_wg,
output_wg,
design_R,
Bottom
]
sim = mp.Simulation(
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
sources=source,
# symmetries=[mp.Mirror(direction=mp.Y)],
default_material=Air,
resolution=resolution,
)
##We can proceed to define our objective function, its corresponding arguments, and the optimization object.
mode = 1
TE0 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3 + 2),
size=mp.Vector3(y=3.5, z = 3),
),
mode,
)
TE_out1 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x= -15),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out2 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x= -12),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out3 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x= -8),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out4 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x= -5),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out5 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x= 0),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out6 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x= 8),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out7 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x= 12),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
ob_list = [TE0, TE_out1, TE_out2, TE_out3, TE_out4, TE_out5, TE_out6, TE_out7]
def J(source, out1, out2, out3, out4, out5, out6, out7):
power1 = npa.abs(out1 / source) ** 2
power2 = npa.abs(out2 / source) ** 2
power3 = npa.abs(out3 / source) ** 2
power4 = npa.abs(out4 / source) ** 2
power5 = npa.abs(out5 / source) ** 2
power6 = npa.abs(out6 / source) ** 2
power7 = npa.abs(out7 / source) ** 2
return npa.mean(power1)
opt = mpa.OptimizationProblem(
simulation=sim,
objective_functions=J,
objective_arguments=ob_list,
design_regions=[design_region],
frequencies=frequencies,
decay_by=1e-5,
)
##We'll define a simple objective function that returns the gradient. We'll plot the new geometry after each iteration.
evaluation_history = []
cur_iter = [0]
def f(v, gradient, cur_beta):
print("Current iteration: {}".format(cur_iter[0] + 1))
f0, dJ_du = opt([v])
plt.figure()
ax = plt.gca()
opt.plot2D(True, output_plane=mp.Volume(size=(np.inf, np.inf, 0), center=(0, 0, 0))) #changed
evaluation_history.append(np.max(np.real(f0)))
#changed
source_coef, out_coef, = opt.get_objective_arguments()
out_profile = np.abs(out_coef / source_coef) ** 2
print("Average through transmission is : ", np.average(out_profile))
cur_iter[0] = cur_iter[0] + 1
return np.real(f0)
algorithm = nlopt.LD_MMA
n = Nx * Ny
##Waveguide
matrix = np.full((Nx, Ny), 0.0)
matrix[:,round((Ny+1)/2 - waveguide_width/2 *resolution) : round((Ny+1)/2 + waveguide_width/2 *resolution)]=1
x = np.reshape(matrix,(n,))
x[Si_mask.flatten()] = 1 # set the edges of waveguides to silicon
x[SiO2_mask.flatten()] = 0 # set the other edges to SiO2
opt.update_design([x])
opt.plot2D(True,output_plane=mp.Volume(center=mp.Vector3(), size=mp.Vector3(np.inf, np.inf, 0)),
field_parameters={'alpha': 0.3})
plt.show()
opt.plot2D(True,output_plane=mp.Volume(center=mp.Vector3(x=0, y=0, z=0), size=mp.Vector3(x = 0, y = 8, z = 8)),
field_parameters={'alpha': 0.3})
plt.title('x = 0')
plt.show()
##lower and upper bounds
lb = np.zeros((Nx * Ny,))
lb[Si_mask.flatten()] = 1
ub = np.ones((Nx * Ny,))
ub[SiO2_mask.flatten()] = 0
cur_beta = 2
beta_scale = 2
num_betas = 7
update_factor = 20
import time
start_time = time.time()
for iters in range(num_betas):
print("current beta: ", cur_beta)
solver = nlopt.opt(algorithm, n)
solver.set_lower_bounds(lb)
solver.set_upper_bounds(ub)
solver.set_max_objective(lambda a, g: f(a, g, cur_beta))
solver.set_maxeval(update_factor)
x[:] = solver.optimize(x)
cur_beta = cur_beta * beta_scale
end_time = time.time()
elapsed_time = end_time - start_time #
print(f"RUNNING TIME：{elapsed_time:.3f} seconds") |
Beta Was this translation helpful? Give feedback.
-
|
Probably you are just computing loss incorrectly. For example, the eigenmode source doesn't exactly launch the desired mode over a nonzero bandwidth, because of modal dispersion (see #2315). Also, your mode monitors are probably too small to properly calculate the mode — they don't include enough of the exponential tails of the mode — since the mode of a ridge waveguide like this is often quite delocalized. Wouldn't it be better to use a mode solver for this? Doing "cutback" experiments in FDTD is a pretty terrible way to do this sort of calculation, since it requires a large 3d calculation when you could be just doing a 2d calculation with a cross section. Since the materials here seem to be lossless, can't you easily obtain a guided mode with zero loss? I don't understand what your optimization objective is here. |
Beta Was this translation helpful? Give feedback.
-
|
Thank you very much for your prompt feedback and suggestions! Based on your suggestions, I also conducted some trials, but did not observe any noticeable effects.1. "I don't understand what your optimization objective is here."Sorry for the ambiguity in my title. My goal is to optimize InP-based devices (e.g., the S-Bend Waveguide, code and optimization logs attached) — not straight waveguides. The S-Bend Waveguide doesn’t seem like a particularly challenging optimization task, as the input/output waveguide centers are spaced only 2 μm apart. While adapting MEEP’s 2D splitter example (04-Splitter.ipynb) to 3D worked well for Si/SiO2, replacing materials with InGaAsP/Air/InP caused non-convergence (According to the logs, the transmission (objective function) only improved from 0.006 to 0.05 over 30+ steps — while the design region’s small size might limit peak transmission, such low values suggest potential model setup errors). The optimization logs show recurring warnings （as shown in log）, and I’m unsure if I’ve overlooked critical information. 2. "For example, the eigenmode source doesn't exactly launch the desired mode over a nonzero bandwidth, because of modal dispersion (see #2315)."I'm not sure if other sources are suitable for subsequent optimization with the Adjoint Solver, I chose the eigenmode source based on inverse design experience (as in (https://github.com/NanoComp/meep/blob/master/python/examples/adjoint_optimization/04-Splitter.ipynb) ). I narrowed the bandwidth of the source (from 1310-1610nm to 1335-1345nm) but observed no improvement. 3. "Wouldn't it be better to use a mode solver for this? Doing 'cutback' experiments in FDTD is a pretty terrible way to do this sort of calculation, since it requires a large 3D calculation when you could be just doing a 2D calculation with a cross-section."My focus is topology optimization of device structures, not modal analysis. While 2D simulations are efficient, they seem unreliable for inverse design outcomes. 3D simulations remain my only option despite computational costs. Do you have any other suggestions? 4. "Also, your mode monitors are probably too small to properly calculate the mode — they don’t include enough of the exponential tails of the mode — since the mode of a ridge waveguide like this is often quite delocalized."Based on your suggestion, I tried different sizes of monitors. I placed different sizes of monitors near the point at 12 microns, and from the results, when the Monitor size was 2, the transmission rate dropped, but when the size was >=3, the transmission rate was almost consistent. This seems to prove that the original code setting the size to 3 is not a problem. In addition, the Monitor size I set in Lumerical FDTD is also (y-span, z-span)=(3,3). Besides, as mentioned before, the ultimate goal is to design a topological structure, so strictly speaking, this straight waveguide is not a Ridge waveguide, InP material serves as the substrate, and the waveguide material is InGaAsP. I am still very puzzled, but I have some guesses：
Thank you again for your expertise. If time permits, could you review my bend waveguide optimization code? Any guidance you can provide would be greatly appreciated. I look forward to your insights! Here is my code for the optimization process of Bend waveguide: import meep as mp
import meep.adjoint as mpa
import autograd.numpy as npa
from autograd import tensor_jacobian_product, grad
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.patches import Circle
import nlopt
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
seed = 240
np.random.seed(seed)
mp.quiet(quietval=True)
InP = mp.Medium(index=3.19)
InP_PQ = mp.Medium(index=3.36)
Air = mp.Medium(index=1)
# As before, we'll define our geometry, filtering, and wavelength parameters.
waveguide_width = 1.5
design_region_width = 5
design_region_height = 5
arm_separation = 1.0
waveguide_length = 2
pml_size = 1.5
resolution = 20
t_InP = 5 #bottom
t_InP_PQ = 0.4 #waveguide
minimum_length = 0.04
eta_e = 0.55
filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)
eta_i = 0.5
eta_d = 1 - eta_e
design_region_resolution = int(5*resolution)
frequencies = 1 / np.linspace(1.4, 1.6, 10)
# We'll also define our simulation domain and set up a broadband source.
Sx = 2 * pml_size + 2 * waveguide_length + design_region_width
Sy = Sx
Sz = 2 * pml_size + t_InP_PQ + 3
cell_size = mp.Vector3(Sx, Sy, Sz)
pml_layers = [mp.PML(pml_size)]
fcen = 0.5 * (1 / 1.4 + 1 / 1.6)
fwidth = 1 / 1.4 - 1 / 1.6
source_center = [-Sx / 2 + pml_size + waveguide_length / 3, 1, 0]
source_size = mp.Vector3(0, 3, 3) #changed
kpoint = mp.Vector3(1, 0, 0)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
mp.EigenModeSource(
src,
eig_band=1,
direction=mp.NO_DIRECTION,
eig_kpoint=kpoint,
size=source_size,
center=source_center,
)
]
# Next, we'll define our design region. This time, however, we'll include a symmetry across the Y plane (normal direction of the symmetry plane points in the Y direction).
Nx = int(design_region_resolution * design_region_width) + 1
Ny = int(design_region_resolution * design_region_height) + 1
design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny), Air, InP_PQ, grid_type="U_MEAN", do_averaging = False) # changed
design_region = mpa.DesignRegion(
design_variables,
volume=mp.Volume(
center=mp.Vector3(),
size=mp.Vector3(design_region_width, design_region_height, t_InP_PQ), #changed
),
)
# We'll define a filtering and interpolation function. We need to include the symmetry requirements in the filter too.
def mapping(x, eta, beta):
# filter
filtered_field = mpa.conic_filter(
x,
filter_radius,
design_region_width,
design_region_height,
design_region_resolution,
)
# projection
projected_field = mpa.tanh_projection(filtered_field, beta, eta)
# interpolate to actual materials
return projected_field.flatten()
# We also need to define a bitmask that describes the boundary conditions of the waveguide and cladding layers.
# Define spatial arrays used to generate bit masks
x_g = np.linspace(-design_region_width / 2, design_region_width / 2, Nx)
y_g = np.linspace(-design_region_height / 2, design_region_height / 2, Ny)
X_g, Y_g = np.meshgrid(x_g, y_g, sparse=True, indexing="ij")
# Define the core mask
left_wg_mask = (X_g == -design_region_width / 2) & (np.abs(Y_g) <= waveguide_width / 2) #INPUT WAVEHUIDE
top_wg_mask = (Y_g == design_region_width / 2) & (np.abs(X_g) <= waveguide_width / 2) #OUTPUT WAVEHUIDE
Si_mask = left_wg_mask | top_wg_mask
# Define the cladding mask
border_mask = (
(X_g == -design_region_width / 2)
| (X_g == design_region_width / 2)
| (Y_g == -design_region_height / 2)
| (Y_g == design_region_height / 2)
)
SiO2_mask = border_mask.copy()
SiO2_mask[Si_mask] = False
S_off = 1
#Finally we can formulate our geometry and simulation object.
input_wg = mp.Block(
center = mp.Vector3(x=-Sx / 4, y=S_off),
material = InP_PQ,
size = mp.Vector3(Sx / 2 + 1, waveguide_width, t_InP_PQ)
)
output_wg = mp.Block(
center = mp.Vector3(x=Sy / 4,y=-S_off),
material = InP_PQ,
size = mp.Vector3(Sx / 2 + 1, waveguide_width, t_InP_PQ)
)
design_R = mp.Block(
center=design_region.center, size=design_region.size, material=design_variables
)
Bottom = mp.Block(
center=mp.Vector3(z = -t_InP_PQ/2-t_InP/2),
material=InP,
size=mp.Vector3(Sx, Sy, t_InP)
)
geometry = [
input_wg,
output_wg,
design_R,
Bottom
]
sim = mp.Simulation(
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
sources=source,
# symmetries=[mp.Mirror(direction=mp.Y)],
default_material=Air,
resolution=resolution,
)
# We can proceed to define our objective function, its corresponding arguments, and the optimization object.
mode = 1
TE0 = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3,y=S_off),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out = mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=mp.Vector3(x=Sx / 2 - pml_size - 2 * waveguide_length / 3,y=-S_off),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out1 = mpa.EigenmodeCoefficient( #for debug
sim,
mp.Volume(
center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3 - 2,y=S_off),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
TE_out2 = mpa.EigenmodeCoefficient( #for debug
sim,
mp.Volume(
center=mp.Vector3(x=Sx / 2 - pml_size - 2 * waveguide_length / 3 + 2,y=-S_off),
size=mp.Vector3(y=3, z = 3),
),
mode,
)
ob_list = [TE0, TE_out, TE_out1, TE_out2]
def J(source, out, out1, out2):
power = npa.abs(out / source) ** 2
return npa.mean(power)
opt = mpa.OptimizationProblem(
simulation=sim,
objective_functions=J,
objective_arguments=ob_list,
design_regions=[design_region],
frequencies=frequencies,
decay_by=1e-5,
)
# We'll define a simple objective function that returns the gradient. We'll plot the new geometry after each iteration.
evaluation_history = []
cur_iter = [0]
def f(v, gradient, cur_beta):
print("Current iteration: {}".format(cur_iter[0] + 1))
f0, dJ_du = opt([mapping(v, eta_i, cur_beta)])
plt.figure()
ax = plt.gca()
opt.plot2D(True, output_plane=mp.Volume(size=(np.inf, np.inf, 0), center=(0, 0, 0))) #changed
circ = Circle((2, 2), minimum_length / 2)
ax.add_patch(circ)
ax.axis("off")
plt.savefig(file_addr + '/run_pic/' + str(cur_iter[0]) + '.png')
# plt.show()
if gradient.size > 0:
gradient[:] = tensor_jacobian_product(mapping, 0)(
v, eta_i, cur_beta, np.sum(dJ_du, axis=1)
)
grad_norm = np.linalg.norm(gradient)
print(f"grad_norm={grad_norm:.8f}")
evaluation_history.append(np.max(np.real(f0)))
#For debug
source_coef, out_coef, TE_out1, TE_out2= opt.get_objective_arguments()
out_profile = np.abs(out_coef / source_coef) ** 2
out_profile1 = np.abs(TE_out1 / source_coef) ** 2
out_profile2 = np.abs(TE_out2 / source_coef) ** 2
print("Average through transmission is : ", np.average(out_profile)," RAW: ",np.average(np.abs(out_coef)), " / ", np.average(np.abs(source_coef)))
print("TE1 transmission is : ", np.average(out_profile1), " RAW: ", np.average(np.abs(TE_out1))) #For debug
print("TE2 transmission is : ", np.average(out_profile2), " RAW: ", np.average(np.abs(TE_out2))) #For debug
#-------
cur_iter[0] = cur_iter[0] + 1
return np.real(f0)
# Finally we'll run the optimizer.
algorithm = nlopt.LD_MMA
n = Nx * Ny # number of parameters
# Initial guess
x = np.ones((n,)) * 0.7
# x = np.random.rand(n)
x[Si_mask.flatten()] = 1 # set the edges of waveguides to silicon
x[SiO2_mask.flatten()] = 0 # set the other edges to SiO2
# lower and upper bounds
lb = np.zeros((Nx * Ny,))
lb[Si_mask.flatten()] = 1
ub = np.ones((Nx * Ny,))
ub[SiO2_mask.flatten()] = 0
cur_beta = 2
beta_scale = 2
num_betas = 7
update_factor = 20 #12
import time
start_time = time.time()
file_addr = "xxxxxxx"
all_x = np.zeros((160,251001))
for iters in range(num_betas):
print("current beta: ", cur_beta)
solver = nlopt.opt(algorithm, n)
solver.set_lower_bounds(lb)
solver.set_upper_bounds(ub)
solver.set_max_objective(lambda a, g: f(a, g, cur_beta))
solver.set_maxeval(update_factor)
x[:] = solver.optimize(x)
cur_beta = cur_beta * beta_scale
end_time = time.time()
elapsed_time = end_time - start_time # Timer
print(f"RUNNING TIME：{elapsed_time:.3f} seconds")Here is the log(p_mp)xx@xxx:~.../python/examples/adjoint_optimization$ PYTHONUNBUFFERED=1 mpirun -np 35 python -m mpi4py Sbend_3D_inP_small.py |
Beta Was this translation helpful? Give feedback.
-
Hello,
I'm attempting to optimize a waveguide structure using the Adjoint solver, where the substrate is a thick InP layer (refractive index = 3.19) and the waveguide consists of InGaAsP (index = 3.36) and air. During optimization, I observed unexpected behavior. To debug, I simplified the designRegion to a straight waveguide and placed multiple mode monitors along its length.
Surprisingly, the measured propagation loss over a 20 μm waveguide segment reaches nearly 1 dB (see Figure 1 for transmittance vs. wavelength at each position). However, the same structure simulated in Lumerical FDTD shows negligible loss. Could you help identify potential errors in my model setup?
mycode.txt
@smartalecH @oskooi @stevengj
Beta Was this translation helpful? Give feedback.
All reactions