Add support for DiffractedPlanewave to EigenmodeCoefficient objective function of adjoint solver
#2054
Conversation
|
Even though we cannot add a unit test in this PR until #2055 is first resolved, here is a draft of the unit test which computes the directional derivative based on an objective function which uses the mode coefficient of a transmitted order of a 2D cylindrical grating with 3D unit cell. The gradient from the adjoint solver is validated using the standard approach of comparing to the brute-force result obtained via a finite difference. Unfortunately, there is a noticeable discrepancy (a factor of ~4) which is likely due to discretization errors that tend to be larger in 3D compared to 2D.
One possibility is to just use a 2D test rather than this 3D test at a larger resolution. import meep as mp
import meep.adjoint as mpa
import math
import numpy as np
from autograd import numpy as npa
from autograd import tensor_jacobian_product
resolution = 25 # pixels/μm
nSi = 3.45
Si = mp.Medium(index=nSi)
nSiO2 = 1.45
SiO2 = mp.Medium(index=nSiO2)
wvl = 0.5 # wavelength
fcen = 1/wvl
dpml = 1.0 # PML thickness
dsub = 3.0 # substrate thickness
dair = 3.0 # air padding
hcyl = 0.5 # cylinder height
rcyl = 0.2 # cylinder radius
sx = 1.1
sy = 0.8
sz = dpml+dsub+hcyl+dair+dpml
cell_size = mp.Vector3(sx,sy,sz)
design_region_size = mp.Vector3(sx,sy)
design_region_resolution = int(2*resolution)
Nx = int(design_region_size.x*design_region_resolution)
Ny = int(design_region_size.y*design_region_resolution)
# ensure reproducible results
rng = np.random.RandomState(59387385)
xcoord = np.linspace(-0.5*design_region_size.x,0.5*design_region_size.x,Nx)
ycoord = np.linspace(-0.5*design_region_size.y,0.5*design_region_size.y,Ny)
xv, yv = np.meshgrid(xcoord,ycoord)
# cylindrical grating
p = np.sqrt(np.square(xv) + np.square(yv)) < rcyl
p = p.flatten(order='F')
# random epsilon perturbation for design region
deps = 1e-5
dp = deps*rng.rand(Nx*Ny)
boundary_layers = [mp.PML(thickness=dpml,direction=mp.Z)]
# periodic boundary conditions
k_point = mp.Vector3()
# plane of incidence: XZ
# polarization: (1) S/TE: Ey or (2) P/TM: Ex
src_cmpt = mp.Ex
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
size=mp.Vector3(sx,sy,0),
center=mp.Vector3(0,0,-0.5*sz+dpml),
component=src_cmpt)]
tran_pt = mp.Vector3(0,0,0.5*sz-dpml)
stop_cond = mp.stop_when_dft_decayed()
substrate = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
material=SiO2)]
def forward_simulation(design_params):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny,0),
mp.air,
Si,
weights=design_params)
matgrid_geometry = [mp.Block(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hcyl),
size=mp.Vector3(design_region_size.x,design_region_size.y,hcyl),
material=matgrid)]
geometry = substrate + matgrid_geometry
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
geometry=geometry,
boundary_layers=boundary_layers,
k_point=k_point)
tran_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)))
sim.run(until_after_sources=stop_cond)
res = sim.get_eigenmode_coefficients(tran_flux,
mp.DiffractedPlanewave([0,0,0],
mp.Vector3(1,0,0),
1,
0))
t_coeffs = res.alpha
tran = abs(t_coeffs[0,0,0])**2
sim.reset_meep()
return tran
def adjoint_solver(design_params):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny,0),
mp.air,
Si,
weights=np.ones((Nx,Ny)))
matgrid_region = mpa.DesignRegion(matgrid,
volume=mp.Volume(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hcyl),
size=mp.Vector3(design_region_size.x,
design_region_size.y,
hcyl)))
matgrid_geometry = [mp.Block(center=matgrid_region.center,
size=matgrid_region.size,
material=matgrid)]
geometry = substrate + matgrid_geometry
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
geometry=geometry,
boundary_layers=boundary_layers,
k_point=k_point)
obj_list = [mpa.EigenmodeCoefficient(sim,
mp.Volume(center=tran_pt,
size=mp.Vector3(sx,sy,0)),
mp.DiffractedPlanewave([0,0,0],
mp.Vector3(1,0,0),
0,
1))]
def J(mode_mon):
return npa.power(npa.abs(mode_mon),2)
opt = mpa.OptimizationProblem(simulation=sim,
objective_functions=J,
objective_arguments=obj_list,
design_regions=[matgrid_region],
frequencies=[fcen])
f, dJ_du = opt([design_params])
sim.reset_meep()
return f, dJ_du
if __name__ == '__main__':
flux_unperturbed = forward_simulation(p)
flux_perturbed = forward_simulation(p+dp)
fd_grad = flux_perturbed-flux_unperturbed
adjsol_obj, adjsol_grad = adjoint_solver(p)
if adjsol_grad.ndim < 2:
adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
adj_scale = (dp[None,:]@adjsol_grad).flatten()
print("obj:, {} (forward simulation), {} (adjoint solver)".format(flux_unperturbed,adjsol_obj))
print("dir-deriv:, {} (finite difference), {} (adjoint solver)".format(fd_grad,adj_scale)) |
I can't think of any intrinsic reason why discretization errors would be larger in 3d than in 2d … except for the fact that it is hard to avoid corner singularities in the fields in 3d, which degrade the discretization accuracy, vs 2d where you don't have these singularities for the Ez polarization. (A way to check this would be to try a 2d simulation with the Hz polarization, in which case you have the same types of singularities.) If corner singularities are a problem, one option would be a smooth test structure in 3d. |
|
I have been working on a unit test in 2D for this feature which can be available now and does not depend on #2055 being fixed first. The test involves computing the transmittance of a binary grating using an oblique incident planewave which is defined using two different but equivalent methods: as (1) an eigenmode number or (2) However, while trying to set up this particular test I discovered that it does not seem to be possible to specify an oblique planewave this way using the import unittest
import meep as mp
import math
import cmath
import numpy as np
class TestDiffractedPlanewave(unittest.TestCase):
@classmethod
def setUp(cls):
cls.resolution = 50 # pixels/μm
cls.dpml = 1.0 # PML thickness
cls.dsub = 3.0 # substrate thickness
cls.dpad = 3.0 # length of padding between grating and PML
cls.gp = 2.6 # grating period
cls.gh = 0.4 # grating height
cls.gdc = 0.5 # grating duty cycle
cls.sx = cls.dpml+cls.dsub+cls.gh+cls.dpad+cls.dpml
cls.sy = cls.gp
cls.cell_size = mp.Vector3(cls.sx,cls.sy,0)
cls.absorber_layers = [mp.PML(thickness=cls.dpml,direction=mp.X)]
cls.wvl = 0.5 # center wavelength
cls.fcen = 1/cls.wvl # center frequency
cls.df = 0.05*cls.fcen # frequency width
cls.ng = 1.5
cls.glass = mp.Medium(index=cls.ng)
cls.eig_parity = mp.ODD_Z
cls.src_pt = mp.Vector3(-0.5*cls.sx+cls.dpml,0,0)
cls.tran_pt = mp.Vector3(0.5*cls.sx-cls.dpml,0,0)
cls.geometry = [mp.Block(material=cls.glass,
size=mp.Vector3(cls.dpml+cls.dsub,mp.inf,mp.inf),
center=mp.Vector3(-0.5*cls.sx+0.5*(cls.dpml+cls.dsub),0,0)),
mp.Block(material=cls.glass,
size=mp.Vector3(cls.gh,cls.gdc*cls.gp,mp.inf),
center=mp.Vector3(-0.5*cls.sx+cls.dpml+cls.dsub+0.5*cls.gh,0,0))]
def run_eigenmode_source(self,m,dp):
ky = m/self.gp
theta = math.asin(ky/(self.fcen*self.ng))
# k (in source medium) with correct length (plane of incidence: XY)
k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(z=1), theta)
print("mode:, {} (m), {:.4f} (theta), ({:.4f},{:.4f},{:.4f}) (k)".format(m,math.degrees(theta),k.x,k.y,k.z))
symmetries = []
if theta == 0:
k = mp.Vector3()
self.eig_parity += mp.EVEN_Y
symmetries = [mp.Mirror(direction=mp.Y)]
if dp:
sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=self.df),
center=self.src_pt,
size=mp.Vector3(0,self.sy,0),
eig_band=mp.DiffractedPlanewave((0,m,0),
mp.Vector3(0,1,0),
1,
0))]
else:
sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=self.df),
center=self.src_pt,
size=mp.Vector3(0,self.sy,0),
direction=mp.AUTOMATIC if theta==0 else mp.NO_DIRECTION,
eig_kpoint=k,
eig_band=1,
eig_parity=self.eig_parity)]
sim = mp.Simulation(resolution=self.resolution,
cell_size=self.cell_size,
boundary_layers=self.absorber_layers,
k_point=mp.Vector3() if dp else k,
geometry=self.geometry,
sources=sources,
symmetries=symmetries)
tran_flux = sim.add_flux(self.fcen,
0,
1,
mp.FluxRegion(center=self.tran_pt,
size=mp.Vector3(0,self.sy,0)))
sim.run(until_after_sources=mp.stop_when_dft_decayed())
tran = mp.get_fluxes(tran_flux)[0]
sim.reset_meep()
return tran
def test_diffracted_planewave(self):
"""Unit test for launching an oblique planewave using an `EigenModeSource`.
Verifies that the transmittance of a binary grating is identical using two
different methods to launch a planewave at oblique incidence.
"""
m = 1
tran_md = self.run_eigenmode_source(m,False)
tran_dp = self.run_eigenmode_source(m,True)
print("tran:, {:.5f} (eigenmode number), {:.5f} (diffracted planewave)".format(tran_md,tran_dp))
self.assertAlmostEqual(tran_md,tran_dp)
if __name__ == '__main__':
unittest.main() |
Actually, this statement is incorrect. Simply flipping the sign of the diffraction order To properly reverse the direction of the planewave (in order to enable the adjoint source), we also need to flip the sign of the component of the wavevector in the non-periodic direction of the cell. This requires some changes to
import meep as mp
import math
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
resolution = 20 # pixels/μm
cell_size = mp.Vector3(14,10,0)
pml_layers = [mp.PML(thickness=2,direction=mp.X)]
fsrc = 1.0 # frequency of planewave (wavelength = 1/fsrc)
n = 1.5 # refractive index of homogeneous material
default_material = mp.Medium(index=n)
m = -7
ky = m/cell_size.y
theta = math.asin(ky/(fsrc*n))
k_point = mp.Vector3(fsrc*n).rotate(mp.Vector3(z=1), theta)
print("stats:, {}, {:.5f}, k=({:.2f},{:.2f})".format(m,math.degrees(theta),k_point.x,k_point.y))
sources = [mp.EigenModeSource(src=mp.ContinuousSource(fsrc),
center=mp.Vector3(),
size=mp.Vector3(y=10),
eig_band=mp.DiffractedPlanewave((0,m,0),
mp.Vector3(0,1,0),
1,
0),
eig_match_freq=True)]
sim = mp.Simulation(cell_size=cell_size,
resolution=resolution,
boundary_layers=pml_layers,
sources=sources,
k_point=mp.Vector3(),
default_material=default_material,
symmetries=[mp.Mirror(mp.Y)] if theta==0 else [])
sim.run(until=100)
nonpml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(10,10,0))
sim.plot2D(fields=mp.Ez,
plot_boundaries_flag=False,
output_plane=nonpml_vol)
if mp.am_master():
plt.title('m = {}, angle = {:.2f}°'.format(m,math.degrees(theta)))
plt.savefig('diffpw_m{}_ang{:.2f}_{}.png'.format(abs(m),
math.degrees(abs(theta)),
'minus' if m < 0 else 'plus'),
bbox_inches='tight') |
|
I think I may have a workaround to needing to use a This is demonstrated in the following example which involves computing the flux into the m=+1 order of a 2D grating with a random design region (black rectangular region in the schematic below). The test involves the usual check of verifying that the gradient of the objective function computed by the adjoint solver matches the brute-force result via finite difference using the directional derivative.
The key is computing the wavevector of the diffracted order manually and then passing this value as the m = 1 # diffraction order in x direction
# wavevector of diffraction order in xz plane
kdiff = mp.Vector3(m/sx,0,(fcen**2-(m/sx)**2)**0.5)As a test, the forward simulation computes the flux in the diffraction order using a 1. forward simulation ( res = sim.get_eigenmode_coefficients(tran_flux,
mp.DiffractedPlanewave([m,0,0],
mp.Vector3(1,0,0),
0,
1))2. adjoint solver (mode number) obj_list = [mpa.EigenmodeCoefficient(sim,
mp.Volume(center=tran_pt,
size=mp.Vector3(sx,sy,0)),
1,
eig_parity=mp.EVEN_Y,
kpoint_func=lambda *not_used: kdiff)]The results demonstrate that these two different approaches produce equivalent results for the objective function (of the unperturbed structure) and the gradient/directional derivative. import meep as mp
import meep.adjoint as mpa
import math
import numpy as np
from autograd import numpy as npa
from autograd import tensor_jacobian_product
resolution = 50 # pixels/μm
nSi = 3.45
Si = mp.Medium(index=nSi)
nSiO2 = 1.45
SiO2 = mp.Medium(index=nSiO2)
theta_d = math.radians(50.0) # deflection angle
wvl = 1.05 # wavelength
fcen = 1/wvl
sx = wvl/math.sin(theta_d) # period in x
sy = 0.5*wvl # period in y
m = 1 # diffraction order in x direction
# wavevector of diffraction order in xz plane
kdiff = mp.Vector3(m/sx,0,(fcen**2-(m/sx)**2)**0.5)
dpml = 1.0 # PML thickness
dsub = 4.0 # substrate thickness
dair = 4.0 # air padding
hblk = 0.5 # block height
xblk = 0.9 # block width in x
yblk = 0.3 # block width in y
sz = dpml+dsub+hblk+dair+dpml
cell_size = mp.Vector3(sx,sy,sz)
design_region_size = mp.Vector3(xblk,yblk,0)
design_region_resolution = int(2*resolution)
Nx = int(design_region_size.x*design_region_resolution)
Ny = int(design_region_size.y*design_region_resolution)
def mapping(w,filter_radius):
filtered_field = mpa.conic_filter(w.reshape((Nx,Ny)),
filter_radius,
design_region_size.x,
design_region_size.y,
design_region_resolution)
return filtered_field.flatten()
# ensure reproducible results
rng = np.random.RandomState(59387385)
# random design region
p = 0.5*rng.rand(Nx*Ny)
fr = 0.1095834 # filter radius
mapped_p = mapping(p,fr)
# random epsilon perturbation for design region
deps = 1e-5
dp = deps*rng.rand(Nx*Ny)
boundary_layers = [mp.PML(thickness=dpml,direction=mp.Z)]
# periodic boundary conditions
k_point = mp.Vector3()
# plane of incidence: XZ
# polarization: (1) S/TE: Ey or (2) P/TM: Ex
src_cmpt = mp.Ex
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
size=mp.Vector3(sx,sy,0),
center=mp.Vector3(0,0,-0.5*sz+dpml),
component=src_cmpt)]
tran_pt = mp.Vector3(0,0,0.5*sz-dpml)
stop_cond = mp.stop_when_dft_decayed()
substrate = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
material=SiO2)]
def forward_simulation(design_params):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
Si,
weights=design_params.reshape(Nx,Ny))
matgrid_geometry = [mp.Block(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hblk),
size=mp.Vector3(design_region_size.x,
design_region_size.y,
hblk),
material=matgrid)]
geometry = substrate + matgrid_geometry
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
geometry=geometry,
boundary_layers=boundary_layers,
k_point=k_point)
tran_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)))
sim.run(until_after_sources=stop_cond)
res = sim.get_eigenmode_coefficients(tran_flux,
mp.DiffractedPlanewave([m,0,0],
mp.Vector3(1,0,0),
0,
1))
t_coeffs = res.alpha
tran = abs(t_coeffs[0,0,0])**2
sim.reset_meep()
return tran
def adjoint_solver(design_params):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny,0),
mp.air,
Si,
weights=np.ones((Nx,Ny)))
matgrid_region = mpa.DesignRegion(matgrid,
volume=mp.Volume(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hblk),
size=mp.Vector3(design_region_size.x,
design_region_size.y,
hblk)))
matgrid_geometry = [mp.Block(center=matgrid_region.center,
size=matgrid_region.size,
material=matgrid)]
geometry = substrate + matgrid_geometry
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
geometry=geometry,
boundary_layers=boundary_layers,
k_point=k_point)
obj_list = [mpa.EigenmodeCoefficient(sim,
mp.Volume(center=tran_pt,
size=mp.Vector3(sx,sy,0)),
1,
eig_parity=mp.EVEN_Y,
kpoint_func=lambda *not_used: kdiff)]
def J(mode_mon):
return npa.power(npa.abs(mode_mon),2)
opt = mpa.OptimizationProblem(simulation=sim,
objective_functions=J,
objective_arguments=obj_list,
design_regions=[matgrid_region],
frequencies=[fcen])
f, dJ_du = opt([design_params])
sim.reset_meep()
return f, dJ_du
flux_unperturbed = forward_simulation(mapped_p)
flux_perturbed = forward_simulation(mapping(p+dp,fr))
fd_grad = flux_perturbed-flux_unperturbed
print("flux:, {}, {}".format(flux_unperturbed,flux_perturbed))
adjsol_obj, adjsol_grad = adjoint_solver(mapped_p)
bp_adjsol_grad = tensor_jacobian_product(mapping,0)(p,fr,adjsol_grad)
if bp_adjsol_grad.ndim < 2:
bp_adjsol_grad = np.expand_dims(bp_adjsol_grad,axis=1)
adj_scale = (dp[None,:]@bp_adjsol_grad).flatten()
print("obj:, {} (forward sim.), {} (adjoint solver))".format(flux_unperturbed,adjsol_obj))
print("dir-deriv:, {} (finite difference), {} (adjoint solver)".format(fd_grad,adj_scale)) |
|
Note that I think the corresponding DiffractedPlanewave adjoint source should (a) flip the simulation |
Looks like it is not currently possible to generate a backward-propagating planewave using this approach. To demonstrate, consider this simple example which involves launching an oblique planewave in a homogeneous medium (adapted from this tutorial): import meep as mp
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
resolution = 20 # pixels/μm
cell_size = mp.Vector3(14,10,0)
pml_layers = [mp.PML(thickness=2,direction=mp.X)]
# rotation angle (in degrees) of planewave
# counter clockwise (CCW) around Z axis, 0° is +X
rot_angle = np.radians(22.3)
fsrc = 1.0 # frequency (wavelength = 1/fsrc)
n = 1.5 # refractive index of homogeneous material
default_material = mp.Medium(index=n)
k_point = -1 * mp.Vector3(fsrc*n).rotate(mp.Vector3(z=1), rot_angle)
sources = [mp.EigenModeSource(src=mp.ContinuousSource(fsrc),
center=mp.Vector3(),
size=mp.Vector3(y=10),
eig_band=mp.DiffractedPlanewave((0,0,0),
mp.Vector3(0,-1,0),
1,
0))]
sim = mp.Simulation(cell_size=cell_size,
resolution=resolution,
boundary_layers=pml_layers,
sources=sources,
k_point=k_point,
default_material=default_material,
symmetries=[mp.Mirror(mp.Y)] if rot_angle==0 else [])
sim.run(until=20)
nonpml_vol = mp.Volume(center=mp.Vector3(),
size=mp.Vector3(10,10,0))
sim.plot2D(fields=mp.Ez,
plot_boundaries_flag=False,
output_plane=nonpml_vol)
if mp.am_master():
plt.savefig('pw.png',bbox_inches='tight')Note that, as suggested in #2069 (comment), we are using the zeroth diffraction ( incorrect The reason this is not working is because the Lines 651 to 653 in 5275f43 change to: else if (k2 > 0)
k[dd - X] = -1 * sqrt(k2);
}This simple modification produces the correct result. correct To fix this bug, the question though is do we want to modify (1) |
|
@oskooi did we decide on a solution here? Seems pretty straightforward to implement |
|
Alternatively it might be just as easy to create a PlanewaveSource from the analytic formula of a planewave. No need to depend on MPB for something this simple. |
|
In fact, for the adjoint solver, I think it's cleaner to separate diffraction orders from eigenmode coefficients. ie have an entirely separate |
Unfortunately, work on this feature is currently blocked by a bug related to oblique modes and |
oskooi
force-pushed
the
diffpw_adjsol
branch
from
9b8bf64
to
4fc6d08
Compare
Codecov Report
❗ Your organization needs to install the Codecov GitHub app to enable full functionality. @@ Coverage Diff @@
## master #2054 +/- ##
==========================================
- Coverage 73.81% 73.72% -0.10%
==========================================
Files 18 18
Lines 5294 5305 +11
==========================================
+ Hits 3908 3911 +3
- Misses 1386 1394 +8
|
|
Hi @smartalecH. I had a chance to revisit this issue recently. I have a workaround which does not use The key to making this work is to specify the The full script is below and demonstrates using a binary grating with a random design region that the directional derivative computed using the adjoint gradient and finite difference agree for S and P polarizations: 1. S polarization 2. P polarization Important: The problem is that this approach does not work for a 2D grating (will post additional results shortly). The broader issue, I think, seems to be related to inaccurate adjoint gradients in 3D calculations in general. See the discussions @mawc2019 and I had putting together the scripts for https://github.com/NanoComp/photonics-opt-testbed/tree/main/Metagrating3D. I am working on 3D tests for #2661 (to be posted shortly as well) which will provide some examples.
"""Validates the adjoint gradient of the diffraction efficiency of a 1D grating.
The 2D test involves computing the gradient of the diffraction efficiency of
the first transmitted order in air of a 1D binary grating on a substrate given
a normally incident planewave from the substrate. The adjoint gradient is
validated using the brute-force finite difference via the directional
derivative.
"""
from typing import Tuple
from autograd import numpy as npa
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import meep as mp
import meep.adjoint as mpa
import numpy as np
RESOLUTION = 100 # pixels/μm
GRATING_PERIOD_UM = 2.5
GRATING_HEIGHT_UM = 0.5
GRATING_DUTY_CYCLE = 0.5
DESIGN_REGION_SIZE = mp.Vector3(
GRATING_HEIGHT_UM, GRATING_DUTY_CYCLE * GRATING_PERIOD_UM, 0
)
DESIGN_REGION_RESOLUTION = int(2 * RESOLUTION)
NX = int(DESIGN_REGION_SIZE.x * DESIGN_REGION_RESOLUTION) + 1
NY = int(DESIGN_REGION_SIZE.y * DESIGN_REGION_RESOLUTION) + 1
def binary_grating_1d(S_pol: bool = True) -> mpa.OptimizationProblem:
"""Sets up the adjoint problem for a 1D binary grating."""
pml_um = 1.0
substrate_um = 3.0
padding_um = 3.0
wavelength_um = 0.5
frequency_center = 1 / wavelength_um
pml_layers = [mp.PML(direction=mp.X, thickness=pml_um)]
n_glass = 1.5
glass = mp.Medium(index=n_glass)
size_x_um = pml_um + substrate_um + GRATING_HEIGHT_UM + padding_um + pml_um
size_y_um = GRATING_PERIOD_UM
cell_size = mp.Vector3(size_x_um, size_y_um, 0)
k_point = mp.Vector3()
eig_parity = mp.ODD_Z if S_pol else mp.EVEN_Z
src_pt = mp.Vector3(-0.5 * size_x_um + pml_um, 0, 0)
sources = [
mp.Source(
mp.GaussianSource(frequency_center, fwidth=0.1 * frequency_center),
component=mp.Ez if S_pol else mp.Hz,
center=src_pt,
size=mp.Vector3(0, size_y_um, 0),
)
]
matgrid = mp.MaterialGrid(
mp.Vector3(NX, NY),
mp.air,
glass,
weights=np.ones((NX, NY)),
)
matgrid_region = mpa.DesignRegion(
matgrid,
volume=mp.Volume(
center=mp.Vector3(
(-0.5 * size_x_um + pml_um + substrate_um +
0.5 * GRATING_HEIGHT_UM),
0,
0
),
size=mp.Vector3(
DESIGN_REGION_SIZE.x,
DESIGN_REGION_SIZE.y,
0
),
),
)
geometry = [
mp.Block(
material=glass,
size=mp.Vector3(pml_um + substrate_um, mp.inf, mp.inf),
center=mp.Vector3(
-0.5 * size_x_um + 0.5 * (pml_um + substrate_um), 0, 0
),
),
mp.Block(
material=matgrid,
size=matgrid_region.size,
center=matgrid_region.center,
),
]
sim = mp.Simulation(
resolution=RESOLUTION,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
sources=sources,
geometry=geometry,
)
tran_pt = mp.Vector3(0.5 * size_x_um - pml_um, 0, 0)
order_y = 1
kdiff = mp.Vector3(
(frequency_center**2 - (order_y / GRATING_PERIOD_UM)**2)**0.5,
order_y / GRATING_PERIOD_UM,
0
)
obj_list = [
mpa.EigenmodeCoefficient(
sim,
mp.Volume(
center=tran_pt,
size=mp.Vector3(0, size_y_um, 0),
),
mode=1,
kpoint_func=lambda *not_used: kdiff,
eig_parity=eig_parity,
),
]
def J(tran_mon):
return npa.abs(tran_mon)**2
opt = mpa.OptimizationProblem(
simulation=sim,
objective_functions=J,
objective_arguments=obj_list,
design_regions=[matgrid_region],
frequencies=[frequency_center],
)
return opt
if __name__ == "__main__":
opt = binary_grating_1d(False)
# ensure reproducible results
rng = np.random.RandomState(9861548)
# random design region
p = 0.9 * rng.rand(NX * NY)
# constant design region
# p = 0.9 * np.ones((NX * NY))
# random perturbation for design region
deps = 1e-5
dp = deps * rng.rand(NX * NY)
f0, dJ_du = opt([p], need_gradient=True)
fig, ax = plt.subplots()
opt.plot2D(init_opt=False, ax=ax)
fig.savefig(
'binary_grating_1d_plot2D.png', dpi=150, bbox_inches='tight'
)
f1, _ = opt([p + dp], need_gradient=False)
if dJ_du.ndim < 2:
dJ_du = np.expand_dims(dJ_du, axis=1)
adj_directional_deriv = (dp[None, :] @ dJ_du).flatten()
fnd_directional_deriv = f1[0] - f0[0]
print(f"grad:, {fnd_directional_deriv} (finite difference), {adj_directional_deriv[0]} (adjoint)") |
Have you tried checking the gradient for a design region that has 1D degrees of freedom, but in a 3D simulation cell? Just to further narrow down the potential bug. Also, did you turn off smoothing in the material grid? |
Fixes #2047.
This is a draft and still needs a unit test.
It seems all that is required to enable this feature is to simply reverse (when necessary) the direction of the wavevector of the
DiffractedPlanewaveobject (itsgproperty) in theplace_adjoint_sourceand__call__functions of theEigenmodeCoefficientclass. No other changes to the adjoint solver are necessary.