Clarity regarding ValueError raised during adjoint optimization

[drinwater]

Hello,

I am a bit of a beginner to the adjoint optimization feature. I am attempting to optimize a 2d grating coupler. However I run into an invalid_argument error from nlopt when running my code. My setup involves an x-z cross-section of a grating coupler with a source at the top. Could you please guide me on what the possible causes for this could be, or if my setup is just entirely wrong?

Thanks

Code

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


# Adapted slightly from symmetric splitter example notebook

# Useful variables
design_width = 10
design_height = 0.22
resolution = 20
cell_size = mp.Vector3(20, 5, 0)

# Materials
Si = mp.Medium(index=3.45)
SiO2 = mp.Medium(index=1.44)

pml_size = 1.0

minimum_length = 0.09
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.5, 1.6, 10)
# print(1/frequencies)

Sx = 20
Sz = 5
cell_size = mp.Vector3(Sx, Sz,0)

pml_layers = [mp.PML(pml_size)]

# Source above design region with vertical direction of propagation
fcen = 1 / 1.56
width = 0.2
fwidth = width * fcen
# source_center = [-Sx / 2 + pml_size + waveguide_length / 3, 0, 0]
source_center = [0, 0.5, 0]
# source_size = mp.Vector3(0, 2, 0)
source_size = mp.Vector3(design_width, 0, 0)
kpoint = mp.Vector3(0, -1, 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,
    )
]

# Design region
Nx = int(design_region_resolution * design_width) + 1
Nz = int(design_region_resolution * design_height) + 1

design_variables = mp.MaterialGrid(mp.Vector3(Nx, Nz, 0), SiO2, Si, grid_type="U_MEAN")
design_region = mpa.DesignRegion(
    design_variables,
    volume=mp.Volume(
        center=mp.Vector3(),
        size=mp.Vector3(design_width, design_height, 0),
    ),
)

# Filtering and interpolaion function with symmetry
def mapping(x, eta, beta):

    # filter
    filtered_field = mpa.conic_filter(
        x,
        filter_radius,
        design_width,
        design_height,
        design_region_resolution,
    )

    # projection
    projected_field = mpa.tanh_projection(filtered_field, beta, eta)

    projected_field = (
        npa.fliplr(projected_field) + projected_field
    ) / 2  # up-down symmetry

    # interpolate to actual materials
    return projected_field.flatten()

# Define spatial arrays used to generate bit masks
x_g = np.linspace(-design_width / 2, design_width / 2, Nx)
z_g = np.linspace(-design_height / 2, design_height / 2, Nz)
X_g, Z_g = np.meshgrid(x_g, z_g, sparse=True, indexing="ij")

# Define the core mask
left_wg_mask = (X_g == -design_width / 2) & (np.abs(Z_g) <= 0.22 / 2)
Si_mask = left_wg_mask

# Define the cladding mask
border_mask = (
    (X_g == -design_width / 2)
    | (X_g == design_width / 2)
    | (Z_g == -design_height / 2)
    | (Z_g == design_height / 2)
)
SiO2_mask = border_mask.copy()
SiO2_mask[Si_mask] = False

# Geometry
geometry = [
    mp.Block(
        center=mp.Vector3(x=-Sx / 4),
        material=Si,
        size=mp.Vector3(Sx / 2 + 1, 0.22,0),
    ),  # left waveguide
    mp.Block(
        center=design_region.center, size=design_region.size, material=design_variables
    ),  # design region
]

sim = mp.Simulation(
    cell_size=cell_size,
    boundary_layers=pml_layers,
    geometry=geometry,
    sources=source,
    symmetries=[mp.Mirror(direction=mp.Y)],
    default_material=SiO2,
    resolution=resolution,
)

mode = 1

# Objective arguments. One at the top of the design region, one in the waveguide.
TE_top = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(x=-8),
        size=mp.Vector3(y=1.5),
    ),
    mode,
)
TE0 = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(
            0, 1, 0
        ),
        size=mp.Vector3(design_width),
    ),
    mode,
    kpoint_func=lambda frequency, mode: mp.Vector3(0, 1, 0)
)
ob_list = [TE0, TE_top]


def J(source, top):
    power = npa.abs(top / source) ** 2
    return power


opt = mpa.OptimizationProblem(
    simulation=sim,
    objective_functions=J,
    objective_arguments=ob_list,
    design_regions=[design_region],
    frequencies=frequencies,
    decay_by=1e-5,
)

plt.figure()
x0 = mapping(
    np.random.rand(
        Nx * Nz,
    ),
    eta_i,
    128,
)
opt.update_design([x0])
opt.plot2D(True)
plt.show()

# %%
evaluation_history = []
cur_iter = [0]
sensitivity = [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(
        False,
        ax=ax,
        plot_sources_flag=False,
        plot_monitors_flag=False,
        plot_boundaries_flag=False,
    )
    circ = Circle((2, 2), minimum_length / 2)
    ax.add_patch(circ)
    ax.axis("off")
    plt.show()

    if gradient.size > 0:
        gradient[:] = tensor_jacobian_product(mapping, 0)(
            v, eta_i, cur_beta, np.sum(dJ_du, axis=1)
        )

    evaluation_history.append(np.max(np.real(f0)))

    sensitivity[0] = dJ_du
    cur_iter[0] = cur_iter[0] + 1

    return np.real(f0)

algorithm = nlopt.LD_MMA
n = Nx * Nz  # number of parameters

# Initial guess
x = np.ones((n,)) * 0.5
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 * Nz,))
lb[Si_mask.flatten()] = 1
ub = np.ones((Nx * Nz,))
ub[SiO2_mask.flatten()] = 0

cur_beta = 4
beta_scale = 2
num_betas = 6
update_factor = 12
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