Example: Planewave Absorption Maximization

In this example, we compute bounds on the maximum absorption of an incident planewave by a structure.

The harmonic time-averaged absorption is

\[f(\mathbf{e}) = \frac{\omega}{2} \int \Im(\epsilon(\mathbf{r})) |\mathbf{e}(\mathbf{r})|^2 d \mathbf{r}\]

Using the fact that \(\Im(\epsilon(\mathbf{r})) |\mathbf{e}(\mathbf{r})|^2\) is only non-zero over the structure, we can rewrite this using the polarization field as

\[f(\mathbf{p}) = \frac{\omega}{2} \Im\left(\chi \frac{\mathbf{p}^\dagger}{\chi^*} \frac{\mathbf{p}}{\chi} \right) = \frac{\omega}{2}\frac{\Im(\chi)}{|\chi|^2} \mathbf{p}^\dagger \mathbf{p}.\]

Therefore, the constant and linear parts of the objective are \(0\), and the quadratic part for the dense form is \(A_0 = -\frac{\omega}{2}\frac{\Im(\chi)}{|\chi|^2} I = \frac{\omega}{2} \Im\left(\frac{1}{\chi} \right) I\). Note the \(-\) sign when defining \(A_0\): this is because the QCQP class writes the QCQP objective as

\[\max_{\mathbf{x}} \quad -\mathbf{x}^\dagger A_0 \mathbf{x} + 2\Re(\mathbf{x^\dagger} s_0) + c_0.\]

import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
import sys, time, os
from dolphindes import photonics, geometry

## wavelength, geometry and materials of the planewave absorption problem ##
wavelength = 1.0
omega = 2*np.pi / wavelength

chi = 3 + 1e-2j

px_per_length = 20 # pixels per length unit. If wavelength = 1.0, then this is pixels per wavelength.
dl = 1 / px_per_length

des_x = 1.5
des_y = 1.5 # size of the design region for the absorbing structure
pmlsep = 1.0
pmlthick = 0.5
Mx = int(des_x / dl)
My = int(des_y / dl)

Npmlsepx = Npmlsepy = int(pmlsep / dl)
Npmlx = Npmly = int(pmlthick / dl)
Nx = Mx + 2*(Npmlsepx + Npmlx)
Ny = My + 2*(Npmlsepy + Npmly)

des_mask = np.zeros((Nx,Ny), dtype=bool)
des_mask[Npmlx+Npmlsepx:-(Npmlx+Npmlsepx) , Npmly+Npmlsepy:-(Npmly+Npmlsepy)] = True
Ndes = int(np.sum(des_mask))

## planewave source
ji = np.zeros((Nx,Ny), dtype=complex)
ji[Npmlx,:] = 2.0 / dl # linesource for unit amplitude planewave traveeling in x direction

plt.matshow(des_mask + np.real(ji)*dl) # visualize where the mask and the source are

## setup geometry
geo = geometry.CartesianFDFDGeometry(
    Nx=Nx, Ny=Ny, Npmlx=Npmlx, Npmly=Npmly, dx=dl, dy=dl
)

## setup the photonic TM FDFD class
abs_problem = photonics.Photonics_TM_FDFD(
    omega=omega, geometry=geo, chi=chi,
    des_mask=des_mask, ji=ji, sparseQCQP=True
)

print(abs_problem)

## obtain incident planewave and plot it
ei = abs_problem.get_ei(ji, update=True)
plt.imshow(np.real(ei), cmap='bwr')
plt.colorbar()

## We start by calculating limits with just the global conservation of power constraints
c0 = 0.0
s0_p = np.zeros(Ndes, dtype=complex)
A0_p = (omega/2) * np.imag(1.0/chi) * sp.eye_array(Ndes)
abs_problem.set_objective(s0=s0_p, A0=A0_p, c0=c0, denseToSparse=True) # using sparse formulation
abs_problem.setup_QCQP(Pdiags='global', verbose=1)
result = abs_problem.bound_QCQP('newton')
print(f'global constraint bound is {result[0]}')

## if desired, we can run generalized constraint descent to tighten the bounds
abs_problem.QCQP.run_gcd()
abs_gcd_bound = abs_problem.QCQP.current_dual

print('Bound for maximum absorption power found is', abs_gcd_bound)
Sx = 0.5 # unit amplitude planewave has Poynting vector magnitude of 0.5
print('Ratio between absorption cross section bound and geometric cross section is', abs_gcd_bound / (Sx*des_y))

Compare with Inverse design

We can perform topology optimization to compare with the bounds. For convenience, dolphindes provides built-in adjoint calculation once the design objective is specified, which can be combined with a non-linear optimization package (like NLopt used here) for easy inverse design.

import nlopt

ndof = int(np.sum(des_mask))
dof = 0.5 * np.ones(ndof) # half slab initialization

opt_data = {'count':0} # dictionary keeping track of optimization progress
def abs_objective(dof, grad):
    # wrapper function around the structure objective in abs_problem to provide printouts of optimization progress
    obj = abs_problem.structure_objective(dof, grad)
    print('at iteration #', opt_data['count'], 'the absorption found is', obj)
    opt_data['count'] += 1
    return obj

# nlopt convergence criteria
maxeval = 1000
ftol_rel = 1e-4

opt = nlopt.opt(nlopt.LD_MMA, ndof)

lb = np.zeros(ndof)
ub = np.ones(ndof)
opt.set_lower_bounds(lb)
opt.set_upper_bounds(ub)
opt.set_ftol_rel(ftol_rel)
opt.set_maxeval(maxeval)

opt.set_max_objective(abs_objective)
opt_dof = opt.optimize(dof)
opt_abs = opt.last_optimum_value()

print('TopOpt found structure with highest absorption of', opt_abs)
print('Ratio between absorption and geometric cross sections is', opt_abs / (Sx*des_y))
print('Compare with the bound ratio:', abs_gcd_bound / (Sx*des_y))