Skip to content

Commit

Permalink
Enable nonzero k_point in adjoint solver (NanoComp#1705)
Browse files Browse the repository at this point in the history 
* fix factor 2

* fix

* fix

* fix

* kpoint

* fix

* Fix the `binary_partition` copy constructor. (NanoComp#1702)

* Fix the `binary_partition` copy constructor.

Prevents issues with older compilers.
\NanoComp#1683 NanoComp#1701

* Update src/structure_dump.cpp

Co-authored-by: Steven G. Johnson <[email protected]>

Co-authored-by: Steven G. Johnson <[email protected]>

* single precision test tol

* single precision test tol

* Fix factor of 2 in adjoint gradients when not using complex fields (NanoComp#1704)

* fix factor 2

* fix

* fix

* fix

* single precision test tol

* single precision test tol

* Update python/adjoint/objective.py

Co-authored-by: Steven G. Johnson <[email protected]>

Co-authored-by: Mo Chen <[email protected]>
Co-authored-by: Steven G. Johnson <[email protected]>

* assertVectorsClose works for scalars as well as vectors (NanoComp#1712)

* assertVectorsClose works for scalars as well as vectors

* rename ApproxComparisonMixin -> ApproxComparisonTestCase, since it is a subclass of TestCase, and use it as such

* add --with-coverage to control usage of Python coverage tests (NanoComp#1713)

* add --with-coverage to control usage of Python coverage tests

* move  --with-coverage  to correct CI line

* flush subnormals on x86 (NanoComp#1709)

* flush subnormals on x86

* less aggressive error threshold in single precision

* increase single-precision tolerance

* update tols in test_array_metadata

* update for assertVectorsClose rename

* lower tolerance in single precision

* Lower tolerance of CW and eigenfrequency solver tests for single precision (NanoComp#1714)

* lower tolerance of CW and eigenfrequency solver tests for single precision

* remove change in default argument from Python wrappers

* fix factor 2

* fix

* kpoint

* fix

* single precision test tol

* single precision test tol

* using real

* using real

* fix rebase

* fix

* fix

Co-authored-by: Mo Chen <[email protected]>
Co-authored-by: Andreas Hoenselaar <[email protected]>
Co-authored-by: Steven G. Johnson <[email protected]>
Co-authored-by: Ardavan Oskooi <[email protected]>
  • Loading branch information
@mochen4 @ahoenselaar
@stevengj @oskooi @mawc2019
5 people authored and mawc2019 committed Nov 3, 2021
1 parent 729ca5d commit 3033f38
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 45 deletions.
2 changes: 1 addition & 1 deletion python/adjoint/objective.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ def _adj_src_scale(self, include_resolution=True):
scale = dV * iomega / adj_src_phase
# compensate for the fact that real fields take the real part of the current,
# which halves the Fourier amplitude at the positive frequency (Re[J] = (J + J*)/2)
if not self.sim.force_complex_fields:
if self.sim.using_real_fields():
scale *= 2
return scale

Expand Down
2 changes: 2 additions & 0 deletions python/adjoint/optimization_problem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -287,6 +287,7 @@ def prepare_adjoint_run(self):
def adjoint_run(self):
# set up adjoint sources and monitors
self.prepare_adjoint_run()
if self.sim.k_point: self.sim.k_point *= -1
for ar in range(len(self.objective_functions)):
# Reset the fields
self.sim.reset_meep()
Expand Down Expand Up @@ -329,6 +330,7 @@ def adjoint_run(self):
self.sim.get_dft_array(dgm, c, f))

# update optimizer's state
if self.sim.k_point: self.sim.k_point *= -1
self.current_state = "ADJ"

def calculate_gradient(self):
Expand Down
18 changes: 9 additions & 9 deletions python/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1928,15 +1928,7 @@ def init_sim(self):
self.fields.require_component(mp.Ez)
self.fields.require_component(mp.Hz)

def use_real(self):
cond1 = self.is_cylindrical and self.m != 0
cond2 = any([s.phase.imag for s in self.symmetries])
cond3 = not self.k_point
cond4 = self.special_kz and self.k_point.x == 0 and self.k_point.y == 0
cond5 = not (cond3 or cond4 or self.k_point == Vector3())
return not (self.force_complex_fields or cond1 or cond2 or cond5)

if use_real(self):
if self.using_real_fields():
self.fields.use_real_fields()
elif verbosity.meep > 0:
print("Meep: using complex fields.")
Expand All @@ -1951,6 +1943,14 @@ def use_real(self):
hook()

self._is_initialized = True

def using_real_fields(self):
cond1 = self.is_cylindrical and self.m != 0
cond2 = any([s.phase.imag for s in self.symmetries])
cond3 = not self.k_point
cond4 = self.special_kz and self.k_point.x == 0 and self.k_point.y == 0
cond5 = not (cond3 or cond4 or self.k_point == Vector3())
return not (self.force_complex_fields or cond1 or cond2 or cond5)

def initialize_field(self, cmpnt, amp_func):
"""
Expand Down
72 changes: 37 additions & 35 deletions python/tests/test_adjoint_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,13 +47,13 @@
fcen = 1/1.55
df = 0.23*fcen
sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),
center=mp.Vector3(-0.5*sxy+dpml,0),
size=mp.Vector3(0,sxy),
center=mp.Vector3(-0.5*sxy+dpml+0.1,0),
size=mp.Vector3(0,sxy-2*dpml),
eig_band=1,
eig_parity=eig_parity)]


def forward_simulation(design_params,mon_type,frequencies=None, use_complex=False):
def forward_simulation(design_params,mon_type,frequencies=None, use_complex=False, k=False):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
silicon,
Expand All @@ -71,13 +71,14 @@ def forward_simulation(design_params,mon_type,frequencies=None, use_complex=Fals
boundary_layers=boundary_layers,
sources=sources,
geometry=geometry,
force_complex_fields=use_complex)
force_complex_fields=use_complex,
k_point = k)
if not frequencies:
frequencies = [fcen]

if mon_type.name == 'EIGENMODE':
mode = sim.add_mode_monitor(frequencies,
mp.ModeRegion(center=mp.Vector3(0.5*sxy-dpml),size=mp.Vector3(0,sxy,0)),
mp.ModeRegion(center=mp.Vector3(0.5*sxy-dpml-0.1),size=mp.Vector3(0,sxy-2*dpml,0)),
yee_grid=True)

elif mon_type.name == 'DFT':
Expand Down Expand Up @@ -108,7 +109,7 @@ def forward_simulation(design_params,mon_type,frequencies=None, use_complex=Fals
return Ez2


def adjoint_solver(design_params, mon_type, frequencies=None, use_complex=False):
def adjoint_solver(design_params, mon_type, frequencies=None, use_complex=False, k=False):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
silicon,
Expand All @@ -129,14 +130,15 @@ def adjoint_solver(design_params, mon_type, frequencies=None, use_complex=False)
boundary_layers=boundary_layers,
sources=sources,
geometry=geometry,
force_complex_fields=use_complex)
force_complex_fields=use_complex,
k_point = k)
if not frequencies:
frequencies = [fcen]

if mon_type.name == 'EIGENMODE':
obj_list = [mpa.EigenmodeCoefficient(sim,
mp.Volume(center=mp.Vector3(0.5*sxy-dpml),
size=mp.Vector3(0,sxy,0)),1)]
mp.Volume(center=mp.Vector3(0.5*sxy-dpml-0.1),
size=mp.Vector3(0,sxy-2*dpml,0)),1)]

def J(mode_mon):
return npa.abs(mode_mon)**2
Expand Down Expand Up @@ -204,31 +206,31 @@ def test_adjoint_solver_DFT_fields(self):

def test_adjoint_solver_eigenmode(self):
print("*** TESTING EIGENMODE ADJOINT FEATURES ***")
## test the single frequency and multi frequency cases
for use_complex in [True, False]:
for frequencies in [[fcen], [1/1.58, fcen, 1/1.53]]:
## compute gradient using adjoint solver
adjsol_obj, adjsol_grad = adjoint_solver(p, MonitorObject.EIGENMODE, frequencies, use_complex)

## compute unperturbed S12
S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies, use_complex)

## compare objective results
print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
self.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)

## compute perturbed S12
S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies, use_complex)

## compare gradients
if adjsol_grad.ndim < 2:
adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
adj_scale = (dp[None,:]@adjsol_grad).flatten()
fd_grad = S12_perturbed-S12_unperturbed
print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
tol = 0.1 if mp.is_single_precision() else 0.03
self.assertClose(adj_scale,fd_grad,epsilon=tol)

## test the single frequency and multi frequency case
for k in [False, mp.Vector3(), mp.Vector3(0.8, 0.6)]:
for use_complex in [True, False]:
for frequencies in [[fcen], [1/1.58, fcen, 1/1.53]]:
## compute gradient using adjoint solver
adjsol_obj, adjsol_grad = adjoint_solver(p, MonitorObject.EIGENMODE, frequencies, use_complex, k)

## compute unperturbed S12
S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies, use_complex, k)

## compare objective results
print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
self.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)

## compute perturbed S12
S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies, use_complex, k)

## compare gradients
if adjsol_grad.ndim < 2:
adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
adj_scale = (dp[None,:]@adjsol_grad).flatten()
fd_grad = S12_perturbed-S12_unperturbed
print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
tol = 0.15 if mp.is_single_precision() else 0.03
self.assertClose(adj_scale,fd_grad,epsilon=tol)

def test_gradient_backpropagation(self):
print("*** TESTING BACKPROP FEATURES ***")
Expand Down Expand Up @@ -267,7 +269,7 @@ def test_gradient_backpropagation(self):
adj_scale = (dp[None,:]@bp_adjsol_grad).flatten()
fd_grad = S12_perturbed-S12_unperturbed
print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
tol = 0.04 if mp.is_single_precision() else 0.03
tol = 0.05 if mp.is_single_precision() else 0.03
self.assertClose(adj_scale,fd_grad,epsilon=tol)


Expand Down

0 comments on commit 3033f38

Please sign in to comment.