Note

Click here to download the full example code or to run this example in your browser via Binder

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

Out:

/opt/conda/lib/python3.11/site-packages/moola/adaptors/dolfin_vector.py:20: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if inner_product is not "custom":

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:

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 0x7fd7638afd70>, 'markers': {'triangle': <dolfin.cpp.mesh.MeshFunctionSizet object at 0x7fd7708dedf0>}}

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()
    gs = fig.add_gridspec(3, 1)
    ax = [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

x0 = np.ones(optimizer.nvar) * 0.5

Optimize!

optimizer.minimize(x0)

Out:

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

global iteration 0
---------------------------------------------
objective = -0.6637293886721712
objective = -0.7275889384562827
objective = -1.7198770360921491
objective = -4.626453892760687
objective = -5.476637222489496
objective = -6.92802534201639
objective = -7.706536582406715
objective = -8.656350620713582
objective = -6.568169480141579
objective = -9.083458152918533
objective = -9.321530500183894
objective = -9.475183556887178
objective = -9.584039511947278
objective = -9.639768654845955
objective = -9.94268761622476
global iteration 1
---------------------------------------------
objective = -9.665770279032362
objective = -9.843405498931912
objective = -10.062486555486593
objective = -10.293206125627993
objective = -10.326892914739474
objective = -10.36101430992344
/builds/gyptis/gyptis.gitlab.io/gyptis/tutorials/optimization/plot_optim.py:190: 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`). Consider using `matplotlib.pyplot.close()`.
  fig = plt.figure(figsize=(4.5, 3))
objective = -10.408196222414823
objective = -10.42077634428563
objective = -10.556359189021462
objective = -7.308428910882467
objective = -10.641483814047184
objective = -9.857372855645474
objective = -10.720820066403544
objective = -10.742125177795131
objective = -10.886217916096026
global iteration 2
---------------------------------------------
objective = -8.78755176531896
objective = -9.6798877243426
objective = -8.4779218936994
objective = -10.614829486209118
objective = -6.740357944287081
objective = -10.559580272836046
objective = -10.929583275677198
objective = -10.948634731448168
objective = -11.096249164844666
objective = -11.12321046184717
objective = -11.200017999548834
objective = -10.411251470969884
objective = -11.263636165520618
objective = -11.284071295736789
objective = -11.290452943706086
global iteration 3
---------------------------------------------
objective = -6.758431801803878
objective = -8.749167396703983
objective = -9.892203698517596
objective = -10.308391072961255
objective = -8.489122815608026
objective = -10.757477719116626
objective = -10.840446636514313
objective = -10.915274070706674
objective = -10.953947307201982
objective = -10.971532627236215
objective = -11.136072331300936
objective = -6.823316564819775
objective = -10.883957493894176
objective = -11.161104653351497
objective = -11.165062473163424
global iteration 4
---------------------------------------------
objective = -9.069046418754203
objective = -6.5661456011369275
objective = -10.260823472452168
objective = -10.54964399736611
objective = -10.637917585332175
objective = -10.671069135622469
objective = -10.687277967780581
objective = -10.698421134862514
objective = -10.706309530816503
objective = -10.711986890914131
objective = -10.716114682337317
objective = -10.719135561437486
objective = -10.72135603862325
objective = -10.722993077004103
objective = -10.724202504048979
global iteration 5
---------------------------------------------
objective = -9.601298915104838
objective = -2.757688906378811
objective = -9.899298097242237
objective = -10.210894417388042
objective = -10.377939519889654
objective = -10.408337822424187
objective = -10.422152795850137
objective = -10.42974481838872
objective = -10.434155289129556
objective = -10.43677122193913
objective = -10.438337897133911
objective = -10.439280939221932
objective = -10.439850193860154
objective = -10.44019437498777
objective = -10.440402671585003
global iteration 6
---------------------------------------------
objective = -9.841494956436081
objective = -2.972016530561061
objective = -9.1435076443665
objective = -10.152274535848406
objective = -10.262866074750235
objective = -10.303337301425179
objective = -10.316729289151807
objective = -10.323514869258942
objective = -10.327091181774254
objective = -10.32900486720834
objective = -10.330036314034547
objective = -10.330594342128892
objective = -10.330896847276524
objective = -10.33106100758912
objective = -10.331150143100745

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

Total running time of the script: ( 7 minutes 15.695 seconds)

Estimated memory usage: 2413 MB

Gallery generated by Sphinx-Gallery