diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 3d6c315c1..8da8db33b 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -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 diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index 6ad79b2f0..31e2522ac 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -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() @@ -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): diff --git a/python/simulation.py b/python/simulation.py index c4eb90563..f7fe7effd 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -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.") @@ -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): """ diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index b8bb7d44c..bd72ccd2f 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -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, @@ -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': @@ -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, @@ -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 @@ -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 ***") @@ -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)