Optimizing a lens#

In this tutorial, we will perform a topology optimization of a lens

We first need to import the Python packages

import matplotlib.colors as mplcolors
import matplotlib.pyplot as plt
import numpy as np

import gyptis as gy

plt.ion()
plt.close("all")

Importantly, we enable automatic differentiation:

gy.use_adjoint(True)

import gyptis.optimize as go
import gyptis.utils as gu

Define EM and geometrical parameters

wavelength = 550  # wavelength of operation (in mm)
Lx, Ly = 4800, 1816.2  # box size
lx, ly = 3000, 250  # lens size
pml_dist = 250
f = 726.5
rtarget = wavelength / 10
yf = -Ly / 2 + pml_dist + ly + f

waist = 1500
position = 0, pml_dist + lx / 2

pmesh = 6

epsilon_min = 1
epsilon_max = 3.48**2

Optimization parameters

filtering_type = "density"
# filtering_type = "sensitivity"
rfilt = ly / 5
maxiter = 15
threshold = (0, 7)

We now define the geometry:

geom = gy.BoxPML(
    dim=2,
    box_size=(Lx, Ly),
    pml_width=(wavelength, wavelength),
)

Now we add a recangular design domain where the permittivity will be optimized

design = geom.add_rectangle(-lx / 2, +pml_dist - Ly / 2, 0, lx, ly)
design, box = geom.fragment(design, geom.box)
# sub = geom.add_rectangle(-Lx / 2, -Ly / 2, 0, Lx, pml_dist)
# sub, box = geom.fragment(sub, box)
target = geom.add_circle(0, yf, 0, rtarget)
target, box = geom.fragment(target, box)

Add physical domains:

geom.add_physical(box, "box")
geom.add_physical(design, "design")
# geom.add_physical(sub, "sub")
geom.add_physical(target, "target")

Out:

7

And set the mesh sizes. A good practice is to have a mesh size that is smaller than the wavelength in the media to resolve the field so size = wavelength / (n*pmesh), with n the refractive index.

geom.set_pml_mesh_size(wavelength / pmesh)
geom.set_size("box", wavelength / pmesh)
geom.set_size("target", wavelength / pmesh)
# geom.set_size("sub", wavelength / (pmesh * np.real(epsilon_max) ** 0.5))
geom.set_size("design", wavelength / (1.2 * pmesh * np.real(epsilon_max) ** 0.5))

geom.build()

Out:

{'mesh': <fenics_adjoint.types.mesh.Mesh object at 0x7fd86b9b2bd0>, 'markers': {'triangle': <dolfin.cpp.mesh.MeshFunctionSizet object at 0x7fd87069a430>}}

We define the incident plane wave. The angle is in radian and theta=0 corresponds to a wave travelling from the bottom.

# gb = gy.GaussianBeam(
#     wavelength=wavelength,
#     angle=0,
#     waist=waist,
#     position=position,
#     Npw=21,
#     dim=2,
#     domain=geom.mesh,
#     degree=2,
# )
gb = gy.PlaneWave(
    wavelength=wavelength,
    angle=0,
    dim=2,
    domain=geom.mesh,
    degree=2,
)

Define the objective function

def objective_function(epsilon_design):
    global simulation

    epsilon = dict(box=1, target=1, design=epsilon_design)

    simulation = gy.Scattering(
        geom,
        epsilon=epsilon,
        source=gb,
        degree=2,
        polarization="TE",
    )

    simulation.solve()

    Hz = simulation.solution["total"]
    nrj = gy.assemble((Hz * Hz.conj).real * simulation.dx("target")) / (
        gy.pi * rtarget**2
    )
    return -nrj

Function to filter and project density

def density_proj_filt(self, density, proj, filt, filtering_type):
    if filtering_type == "density":
        density_f = self.filter.apply(density) if filt else density
    else:
        density_f = density
    density_fp = (
        go.projection(density_f, beta=2**self.proj_level) if proj else density_f
    )
    fs_plot = gy.dolfin.FunctionSpace(self.submesh, "DG", 0)
    return gu.project_iterative(density_fp, fs_plot)

Define callback function

jopt = 0


def callback(self):
    global jopt
    proj = self.proj_level is not None
    filt = self.filter != 0
    density_fp = density_proj_filt(self, self.density, proj, filt, filtering_type)

    Hz = simulation.solution["total"]
    projspace = optimizer.fs_ctrl
    # projspace = simulation.real_function_space
    fieldplot = gu.project_iterative((Hz * Hz.conj).real, projspace)
    fig = plt.figure(figsize=(4.5, 3))
    plt.clf()
    ax = []
    gs = fig.add_gridspec(3, 1)
    ax.append(fig.add_subplot(gs[0:2, :]))

    gy.plot(
        fieldplot,
        ax=ax[0],
        cmap="inferno",
        edgecolors="face",
        vmin=0,
        vmax=12,
        # norm=mplcolors.LogNorm(vmin=1e-2, vmax=None),
    )
    geom_lines = geom.plot_subdomains(color="white", ax=ax[0])
    gy.dolfin.plot(density_fp, vmin=0, vmax=1, cmap="Greys", alpha=0.7)
    plt.xlabel(r"$x$ (μm)")
    plt.ylabel(r"$y$ (μm)")
    plt.axis("off")
    # plt.tight_layout()
    ax.append(fig.add_subplot(gs[2, :]))
    plt.sca(ax[1])

    dplot, cb2 = gy.plot(density_fp, ax=ax[1], vmin=0, vmax=1, cmap="Reds", alpha=1)
    plt.xlim(-lx / 2, lx / 2)
    y0 = -Ly / 2 + pml_dist
    plt.ylim(y0, y0 + ly)
    plt.axis("off")
    plt.suptitle(f"iteration {jopt}, objective = {-self.objective:.5f}")
    plt.tight_layout()
    gy.pause(0.1)

    jopt += 1
    return self.objective

Initialize the optimizer

optimizer = go.TopologyOptimizer(
    objective_function,
    geom,
    eps_bounds=(epsilon_min, epsilon_max),
    rfilt=rfilt,
    filtering_type=filtering_type,
    callback=callback,
    verbose=True,
    ftol_rel=1e-6,
    xtol_rel=1e-12,
    maxiter=maxiter,
    threshold=threshold,
)

Initial value

Optimize!

optimizer.minimize(x0)
  • iteration 0, objective = 0.66373
  • iteration 1, objective = 0.72759
  • iteration 2, objective = 1.71988
  • iteration 3, objective = 4.62645
  • iteration 4, objective = 5.47664
  • iteration 5, objective = 6.92803
  • iteration 6, objective = 7.70654
  • iteration 7, objective = 8.65635
  • iteration 8, objective = 6.56817
  • iteration 9, objective = 9.08346
  • iteration 10, objective = 9.32153
  • iteration 11, objective = 9.47518
  • iteration 12, objective = 9.58404
  • iteration 13, objective = 9.63977
  • iteration 14, objective = 9.94269
  • iteration 15, objective = 9.66577
  • iteration 16, objective = 9.84341
  • iteration 17, objective = 10.06249
  • iteration 18, objective = 10.29321
  • iteration 19, objective = 10.32689
  • iteration 20, objective = 10.36101
  • iteration 21, objective = 10.40820
  • iteration 22, objective = 10.42078
  • iteration 23, objective = 10.55636
  • iteration 24, objective = 7.30843
  • iteration 25, objective = 10.64148
  • iteration 26, objective = 9.85737
  • iteration 27, objective = 10.72082
  • iteration 28, objective = 10.74213
  • iteration 29, objective = 10.88622
  • iteration 30, objective = 8.78755
  • iteration 31, objective = 9.67989
  • iteration 32, objective = 8.47792
  • iteration 33, objective = 10.61483
  • iteration 34, objective = 6.74036
  • iteration 35, objective = 10.55958
  • iteration 36, objective = 10.92958
  • iteration 37, objective = 10.94863
  • iteration 38, objective = 11.09625
  • iteration 39, objective = 11.12321
  • iteration 40, objective = 11.20002
  • iteration 41, objective = 10.41125
  • iteration 42, objective = 11.26364
  • iteration 43, objective = 11.28407
  • iteration 44, objective = 11.29045
  • iteration 45, objective = 6.75843
  • iteration 46, objective = 8.74917
  • iteration 47, objective = 9.89220
  • iteration 48, objective = 10.30839
  • iteration 49, objective = 8.48912
  • iteration 50, objective = 10.75748
  • iteration 51, objective = 10.84045
  • iteration 52, objective = 10.91527
  • iteration 53, objective = 10.95395
  • iteration 54, objective = 10.97153
  • iteration 55, objective = 11.13607
  • iteration 56, objective = 6.82332
  • iteration 57, objective = 10.88396
  • iteration 58, objective = 11.16110
  • iteration 59, objective = 11.16506
  • iteration 60, objective = 9.06905
  • iteration 61, objective = 6.56615
  • iteration 62, objective = 10.26082
  • iteration 63, objective = 10.54964
  • iteration 64, objective = 10.63792
  • iteration 65, objective = 10.67107
  • iteration 66, objective = 10.68728
  • iteration 67, objective = 10.69842
  • iteration 68, objective = 10.70631
  • iteration 69, objective = 10.71199
  • iteration 70, objective = 10.71611
  • iteration 71, objective = 10.71914
  • iteration 72, objective = 10.72136
  • iteration 73, objective = 10.72299
  • iteration 74, objective = 10.72420
  • iteration 75, objective = 9.60130
  • iteration 76, objective = 2.75769
  • iteration 77, objective = 9.89930
  • iteration 78, objective = 10.21089
  • iteration 79, objective = 10.37794
  • iteration 80, objective = 10.40834
  • iteration 81, objective = 10.42215
  • iteration 82, objective = 10.42974
  • iteration 83, objective = 10.43416
  • iteration 84, objective = 10.43677
  • iteration 85, objective = 10.43834
  • iteration 86, objective = 10.43928
  • iteration 87, objective = 10.43985
  • iteration 88, objective = 10.44019
  • iteration 89, objective = 10.44040
  • iteration 90, objective = 9.84149
  • iteration 91, objective = 2.97202
  • iteration 92, objective = 9.14351
  • iteration 93, objective = 10.15227
  • iteration 94, objective = 10.26287
  • iteration 95, objective = 10.30334
  • iteration 96, objective = 10.31673
  • iteration 97, objective = 10.32351
  • iteration 98, objective = 10.32709
  • iteration 99, objective = 10.32900
  • iteration 100, objective = 10.33004
  • iteration 101, objective = 10.33059
  • iteration 102, objective = 10.33090
  • iteration 103, objective = 10.33106
  • iteration 104, objective = 10.33115

Out:

#################################################
Topology optimization with 3708 variables
#################################################

global iteration 0
---------------------------------------------
objective = -0.6637293886726063
objective = -0.7275889384561683
objective = -1.7198770360951376
objective = -4.626453892741269
objective = -5.4766372224679705
objective = -6.9280253420767615
objective = -7.70653658245667
objective = -8.656350620751398
objective = -6.568169480228674
objective = -9.083458152957547
objective = -9.3215305002255
objective = -9.47518355694992
objective = -9.584039511999082
objective = -9.639768654902696
objective = -9.942687616276055
global iteration 1
---------------------------------------------
objective = -9.665770279093655
objective = -9.843405498995548
objective = -10.062486555573262
objective = -10.293206125677099
objective = -10.326892914739155
objective = -10.361014309995955
/builds/gyptis/gyptis.gitlab.io/gyptis/tutorials/optimization/plot_optim.py:189: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig = plt.figure(figsize=(4.5, 3))
objective = -10.408196222452416
objective = -10.420776344321785
objective = -10.556359189053612
objective = -7.308428907767103
objective = -10.641483813865738
objective = -9.857372854769011
objective = -10.72082006642754
objective = -10.742125177800494
objective = -10.88621791597908
global iteration 2
---------------------------------------------
objective = -8.78755176519255
objective = -9.679887724218244
objective = -8.477921893501769
objective = -10.614829486056989
objective = -6.740357945075127
objective = -10.559580272477183
objective = -10.929583275555514
objective = -10.948634731267942
objective = -11.096249164765206
objective = -11.123210461765675
objective = -11.200017999568752
objective = -10.411251472427505
objective = -11.26363616546792
objective = -11.284071295680457
objective = -11.290452943654058
global iteration 3
---------------------------------------------
objective = -6.7584318017476495
objective = -8.749167396823127
objective = -9.892203698434399
objective = -10.308391072923174
objective = -8.489122815177396
objective = -10.757477719086324
objective = -10.840446636474798
objective = -10.915274070680553
objective = -10.953947307190768
objective = -10.971532627221391
objective = -11.13607233125303
objective = -6.823316560664409
objective = -10.88395749344119
objective = -11.161104653280214
objective = -11.165062473299772
global iteration 4
---------------------------------------------
objective = -9.06904641860111
objective = -6.56614560093964
objective = -10.260823472261967
objective = -10.549643997181532
objective = -10.637917585162095
objective = -10.671069135457566
objective = -10.687277967616623
objective = -10.698421134699531
objective = -10.706309530653034
objective = -10.711986890750422
objective = -10.716114682173155
objective = -10.719135561273173
objective = -10.721356038460224
objective = -10.72299307684145
objective = -10.724202503885683
global iteration 5
---------------------------------------------
objective = -9.601298914962397
objective = -2.7576889065724517
objective = -9.89929809704211
objective = -10.210894417104736
objective = -10.377939519685137
objective = -10.408337822253298
objective = -10.422152795686703
objective = -10.42974481822942
objective = -10.43415528897249
objective = -10.436771221783529
objective = -10.43833789697882
objective = -10.439280939067087
objective = -10.439850193705077
objective = -10.440194374833114
objective = -10.440402671431436
global iteration 6
---------------------------------------------
objective = -9.841494956248027
objective = -2.9720165314273577
objective = -9.143507644127821
objective = -10.15227453568833
objective = -10.262866074632358
objective = -10.3033373013252
objective = -10.316729289054326
objective = -10.323514869162336
objective = -10.327091181678286
objective = -10.329004867112115
objective = -10.330036313938665
objective = -10.33059434203291
objective = -10.330896847180599
objective = -10.331061007493277
objective = -10.331150143004697

(array([0.3579, 0.4559, 0.3385, ..., 0.0158, 0.493 , 0.005 ]), -10.331150143004697)

Total running time of the script: ( 11 minutes 38.150 seconds)

Estimated memory usage: 2483 MB

Gallery generated by Sphinx-Gallery