Fields blowing up 3D topology optimization

[AhmedHBayoumi]

Hi,

While expanding the 2D splitter example to 3D, the fields started to blow up. I know this issue is mentioned in the Q&A here, but I think none of these cases apply to the splitter example. However, I started to tweak some parameters (some of which are just out of despair) to see if it has an effect without succeeding so far. Here are the things I tried:

1- Thicker PML (2 µm).
2- Using Absorber instead of PML.
3- setting ''' eps_averaging ''' to ''' False '''
4- Increasing the resolution up to 50.
5- Setting '''do_averaging ''' to '''False'''
6- Decreasing Courant factor

This is how I built my 3D model:
1- cell_size is now 3D
2- source and monitors are now 2D planes, with source having a k vector in x direction.
3- MaterialGrid is 2D and DesignRegion is 3D.
4- Waveguide blocks are now 3D.

This is how the geometry looks like:

First iteration gives me 400% transmission at each port, and here's the gradient:

For this geometry:

Which looks suspicious to me, specially that sometimes I get some 'bumps' in the gradients as in here:

Any leads on what I should change or test would be greatly appreciated.

Thanks

Ahmed

[smartalecH]

In your material grid, set do_averaging=False

[smartalecH]

Oh I see you've tried that. Hmmm

[smartalecH]

Can you post a gist of your script? I suspect there's an issue with the setup not obvious from the description.

[AhmedHBayoumi]

Can you post a gist of your script? I suspect there's an issue with the setup not obvious from the description.

Sure, here you are:

Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.44)
waveguide_width = 0.5
design_region_width = 3
design_region_height = 3
arm_separation = 2
waveguide_length = 0.5
t_Si = 0.21
pml_size = 3
resolution = 15
t_oxide = 0.5
minimum_length = 0.08
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)
Sx = 2 * pml_size + 2 * waveguide_length + design_region_width +4
Sy = 2 * pml_size + design_region_height + 4
Sz = 2* pml_size + 5*t_oxide + t_Si + 2
cell_size = mp.Vector3(Sx, Sy, Sz)
pml_layers = [mp.PML(pml_size)]

fcen = 1 / 1.56
width = 0.2
fwidth = width * fcen
source_center = [-design_region_width/2 - 1  , 0, 0]
source_size = mp.Vector3(0, 1.8, 1.8)
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,
    )
]


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), SiO2, Si, grid_type="U_MEAN", do_averaging = False)
design_region = mpa.DesignRegion(
    design_variables,
    volume=mp.Volume(
        center=mp.Vector3(),
        size=mp.Vector3(design_region_width, design_region_height, t_Si),
    ),
)


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)

    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_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)
top_right_wg_mask = (X_g == design_region_width / 2) & (
    np.abs(Y_g + arm_separation / 2) <= waveguide_width / 2
)
bottom_right_wg_mask = (X_g == design_region_width / 2) & (
    np.abs(Y_g - arm_separation / 2) <= waveguide_width / 2
)
Si_mask = left_wg_mask | top_right_wg_mask | bottom_right_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

geometry = [
    mp.Block(
        center=mp.Vector3(x=-Sx / 4),
        material=Si,
        size=mp.Vector3(Sx / 2 + 1, waveguide_width, t_Si),
    ),  # left waveguide
    mp.Block(
        center=mp.Vector3(x=Sx / 4, y=arm_separation / 2),
        material=Si,
        size=mp.Vector3(Sx / 2 + 1, waveguide_width, t_Si),
    ),  # top right waveguide
    mp.Block(
        center=mp.Vector3(x=Sx / 4, y=-arm_separation / 2),
        material=Si,
        size=mp.Vector3(Sx / 2 + 1, waveguide_width, t_Si),
    ),  # bottom right waveguide
    mp.Block(
        center=design_region.center, size=design_region.size, material=design_variables
    ),
]

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,
    eps_averaging = True,
)

mode = 1

TE0 = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(x=-design_region_width/2 -0.5),
        size=mp.Vector3(y=1.8 , z= 1.8),
    ),
    mode,
)
TE_top = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(
            design_region_width/2 + 1, arm_separation/2, 0
        ),
        size=mp.Vector3(y=1.8, z = 1.8),
    ),
    mode,
)
TE_bottom = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(
            design_region_width/2 + 1, -arm_separation / 2, 0
        ),
        size=mp.Vector3(y=1.8, z = 1.8),
    ),
    mode,
)
ob_list = [TE0, TE_top, TE_bottom]


def J(source, top, bottom):
    power = npa.abs(top / source) ** 2 + npa.abs(bottom / 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,
)


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)))
    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)))
    source_coef, top_coef, bottom_ceof = opt.get_objective_arguments()
    top_profile = np.abs(top_coef / source_coef) ** 2
    bottom_profile = np.abs(bottom_ceof / source_coef) ** 2

    print("Average through transmission is : ", np.average(top_profile))
    print("Average cross transmission is : ", np.average(bottom_profile))
       
    plt.figure()
    backprop_gradient = tensor_jacobian_product(mapping,0)(v, eta_i,cur_beta,dJ_du[:,0])
    plt.imshow(np.rot90(backprop_gradient.reshape(Nx,Ny)))
    plt.colorbar()    
    plt.show()
    cur_iter[0] = cur_iter[0] + 1

    return np.real(f0)


algorithm = nlopt.LD_MMA
n = Nx * Ny  # 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 * Ny,))
lb[Si_mask.flatten()] = 1
ub = np.ones((Nx * Ny,))
ub[SiO2_mask.flatten()] = 0

cur_beta = 8
beta_scale = 4
num_betas = 10
update_factor = 30
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

[smartalecH]

Try running it without any symmetries.

[AhmedHBayoumi]

The first iteration seems to be going in the right direction and the values do make sense:

Does that mean that the symmetry could've resulted in a divergence in the fields? How's that possible?

I will update tomorrow when more iterations are done if both FOM and gradient evolution are logical, thanks a lot for your help.

[smartalecH]

So there are two nuances here.

First, there's a known bug when computing adjoint fields with symmetries: #2145. Specifically, only half the fields get computed. (You can manually get around this by mirroring your gradient array until this is fixed).

But that's not what's going on here. Chances are you are specifying something (either with the geometry or the fields/sources) that no longer obeys the symmetry, so the resulting simulation is "undefined". The phase of the symmetry is particularly important. I didn't look too closely to see where that is in your code, but it's a common mistake.

[AhmedHBayoumi]

So indeed taking off the symmetry solves it:

Even though I can't directly figure what was not symmetric in nature in the simulation, but I guess that could be the discretization of the grid as it was showing warning of non integer number of pixels.

Thanks a lot, Alec!