Animate2D for OptimizationProblem

[thomasdorch]

The purpose of this PR is to allow plotting of topology optimizations (#2185).

So far, I've modified the Animate2D class to support changing epsilon values in addition to field values.

[thomasdorch]

Below is a modified version of the Filtering and Thresholding Design Parameters and Broadband Objective Functions tutorial that I'm using to test this.

Currently, I'm just adding the Animate2D object to the function, f, so that it gets called on every iteration.
This now works by passing the Animate2D object into OptimizationProblem within a step function with the added step_funcs argument.

import meep as mp
import meep.adjoint as mpa
import nlopt
import numpy as np
from autograd import numpy as npa
from autograd import tensor_jacobian_product
from matplotlib import pyplot as plt
from matplotlib.patches import Circle
from meep import Animate2D

mp.Verbosity(0)

Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.44)

# Let's choose our basic geometry parameters, along with the frequency range to simulate. If our frequency points are too close together than the resulting adjoint simulation may take longer to simulate each iteration.
waveguide_width = 0.5
design_region_width = 2.5
design_region_height = 2.5

waveguide_length = 0.5
pml_size = 1.0
resolution = 20
frequencies = 1 / np.linspace(1.5, 1.6, 10)

minimum_length = 0.09  # minimum length scale (microns)
eta_i = (
    0.5  # blueprint (or intermediate) design field thresholding point (between 0 and 1)
)
eta_e = 0.55  # erosion design field thresholding point (between 0 and 1)
eta_d = 1 - eta_e  # dilation design field thresholding point (between 0 and 1)
filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)
print(filter_radius)
design_region_resolution = int(1 * resolution)


# As before, we'll generate our simulation cell with the specified sources, boundary layers, geometry, etc.
Sx = 2 * pml_size + 2 * waveguide_length + design_region_width
Sy = 2 * pml_size + design_region_height + 0.5
cell_size = mp.Vector3(Sx, Sy)

pml_layers = [mp.PML(pml_size)]

fcen = 1 / 1.55
width = 0.2
fwidth = width * fcen
source_center = [-Sx / 2 + pml_size + waveguide_length / 3, 0, 0]
source_size = mp.Vector3(0, Sy, 0)
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)
Ny = int(design_region_resolution * design_region_height)

design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny), SiO2, Si, grid_type="U_MEAN")
design_region = mpa.DesignRegion(
    design_variables,
    volume=mp.Volume(
        center=mp.Vector3(),
        size=mp.Vector3(design_region_width, design_region_height, 0),
    ),
)

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")

left_wg_mask = (X_g == -design_region_width / 2) & (np.abs(Y_g) <= waveguide_width / 2)
top_wg_mask = (Y_g == design_region_width / 2) & (np.abs(X_g) <= waveguide_width / 2)
Si_mask = left_wg_mask | top_wg_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

def mapping(x, eta, beta):
    x = npa.where(Si_mask.flatten(), 1, npa.where(SiO2_mask.flatten(), 0, x))
    # 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()

geometry = [
    mp.Block(
        center=mp.Vector3(x=-Sx / 4), material=Si, size=mp.Vector3(Sx / 2, 0.5, 0)
    ),  # horizontal waveguide
    mp.Block(
        center=mp.Vector3(y=Sy / 4), material=Si, size=mp.Vector3(0.5, Sy / 2, 0)
    ),  # vertical waveguide
    mp.Block(
        center=design_region.center, size=design_region.size, material=design_variables
    ),  # design region
    mp.Block(
        center=design_region.center,
        size=design_region.size,
        material=design_variables,
        e1=mp.Vector3(x=-1).rotate(mp.Vector3(z=1), np.pi / 2),
        e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), np.pi / 2),
    ),
]

sim = mp.Simulation(
    cell_size=cell_size,
    boundary_layers=pml_layers,
    geometry=geometry,
    sources=source,
    default_material=SiO2,
    resolution=resolution,
)

mode = 1

TE0 = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(x=-Sx / 2 + pml_size + 2 * waveguide_length / 3),
        size=mp.Vector3(y=Sy),
    ),
    mode,
)
TE_top = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mp.Vector3(0, Sx / 2 - pml_size - 2 * waveguide_length / 3, 0),
        size=mp.Vector3(x=Sx),
    ),
    mode,
)

ob_list = [TE0, TE_top]

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

# Create the animation and specify update_epsilon=True
anim = Animate2D(
    fields=mp.Ez,
    realtime=True,
    update_epsilon=True,
    field_parameters={"alpha": 0.6, "cmap": "RdBu", "interpolation": "none"},
    eps_parameters={'contour': False, 'alpha': 1, 'frequency': 1/1.55},
)

opt = mpa.OptimizationProblem(
    simulation=sim,
    objective_functions=J,
    objective_arguments=ob_list,
    design_regions=[design_region],
    frequencies=frequencies,
    step_funcs=[mp.at_every(25, anim)]
)

evaluation_history = []
cur_iter = [0]

def f(v, gradient, cur_beta):
    print(f"Current iteration: {cur_iter[0] + 1}")

    f0, dJ_du = opt([mapping(v, eta_i, cur_beta)])  # compute objective and gradient

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

    evaluation_history.append(np.real(f0))

    cur_iter[0] = cur_iter[0] + 1

    return np.real(f0)


# We'll now run our optimizer in loop. The loop will increase beta and reset the optimizer, which is important since the cost function changes.

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

for iters in range(num_betas):
    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

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 24.00000% with 38 lines in your changes missing coverage. Please review.

Project coverage is 72.93%. Comparing base (2aa9164) to head (36cf445).
Report is 210 commits behind head on master.

Files with missing lines	Patch %	Lines
python/visualization.py	22.44%	38 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2186      +/-   ##
==========================================
- Coverage   73.09%   72.93%   -0.17%     
==========================================
  Files          17       17              
  Lines        4906     4925      +19     
==========================================
+ Hits         3586     3592       +6     
- Misses       1320     1333      +13     

Files with missing lines	Coverage Δ
python/adjoint/optimization_problem.py	57.06% <100.00%> (+0.22%)	⬆️
python/visualization.py	41.49% <22.44%> (-1.16%)	⬇️

... and 2 files with indirect coverage changes

[stevengj]

Is this ready for review? If so, change the status from "draft".

[thomasdorch]

I wanted to get it a little better integrated into OptimizationProblem before marking it for review. With the step_funcs addition I think it's ready now though.

[thomasdorch]

These tests aren't failing on my end..

[Yaraslaut]

These tests aren't failing on my end..

I have same errors in my PR #2188. My guess that this is introduced by some library update, and valgrind shows me that error comes from mpb.cpp

[thomasdorch]

These tests aren't failing on my end..

I have same errors in my PR #2188. My guess that there are introduce my some library update, and valgrind shows me that error comes from mpb.cpp

Looks like the CI is failing for the master branch as well..

They started 3 days ago, which coincides with the latest commit in MPB, any chance that's related? I can't get the master branch of MPB to pass its own tests, but the v1.11.1 tag works fine.
@Yaraslaut what version of MPB are you on?

[Yaraslaut]

These tests aren't failing on my end..

I have same errors in my PR #2188. My guess that there are introduce my some library update, and valgrind shows me that error comes from mpb.cpp

Looks like the CI is failing for the master branch as well..

They started 3 days ago, which coincides with the latest commit in MPB, any chance that's related? I can't get the master branch of MPB to pass its own tests, but the v1.11.1 tag works fine. @Yaraslaut what version of MPB are you on?

I was using HEAD of the repo to link with meep, I was also getting errors in make check but it looks like this errors come from guile and not from library itself. Ok, I will run tests on another mpb commit then and see what happens.

Update: @thomasdorch you are right, I checked with another mpb version and everything works fine, except test_mpb.py

[oskooi]

LGTM.

Would you post an example just to see what the output actually looks like?

[thomasdorch]

These aren't spaced evenly in time but they illustrate a few snapshots as the optimization progresses:

With fields=mp.Ez:

With fields=None:



... and so on

[smartalecH]

Very cool! Wish we did this earlier...

[stevengj]

LGTM. Can you add documentation, e.g. at demonstrate this in a tutorial?

[thomasdorch]

I am working on getting the animation working in Jupyter notebooks so once that works I can just replace some of the current plot2D calls in the adjoint example notebooks with the animation. How does that sound?