From 0cb3213f1f63847b710d734bb1edb2179ad3b87f Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Fri, 5 Aug 2022 12:16:15 -0500 Subject: [PATCH 01/10] Implemenation mostly typed up, but not yet working due to a bug. --- python/adjoint/objective.py | 121 +++++++++++++++++++++++++++- python/tests/test_adjoint_solver.py | 41 +++++++++- 2 files changed, 158 insertions(+), 4 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 07aeea618..1ae128327 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -10,6 +10,8 @@ from .filter_source import FilteredSource Grid = namedtuple("Grid", ["x", "y", "z", "w"]) +E_CPTS = [mp.Ex, mp.Ey, mp.Ez] +H_CPTS = [mp.Hx, mp.Hy, mp.Hz] class ObjectiveQuantity(abc.ABC): @@ -268,12 +270,14 @@ def __call__(self): class FourierFields(ObjectiveQuantity): - def __init__(self, sim, volume, component, yee_grid=False, decimation_factor=0): + def __init__(self, sim, volume, component, yee_grid=False, decimation_factor=0,scale = 1,conj = False): super().__init__(sim) self.volume = sim._fit_volume_to_simulation(volume) self.component = component self.yee_grid = yee_grid self.decimation_factor = decimation_factor + self.scale = scale + self.conj = conj def register_monitors(self, frequencies): self._frequencies = np.asarray(frequencies) @@ -316,10 +320,13 @@ def place_adjoint_source(self, dJ): self.all_fouriersrcdata = self._monitor.swigobj.fourier_sourcedata( self.volume.swigobj, self.component, self.sim.fields, dJ ) - + for fourier_data in self.all_fouriersrcdata: amp_arr = np.array(fourier_data.amp_arr).reshape(-1, self.num_freq) - scale = amp_arr * self._adj_src_scale(include_resolution=False) + scale = self.scale * amp_arr * self._adj_src_scale(include_resolution=False) + + if self.conj: + scale = np.conjugate(scale) if self.num_freq == 1: sources += [ @@ -483,3 +490,111 @@ def __call__(self): self.ldos_scale = self.sim.ldos_scale self.ldos_Jdata = self.sim.ldos_Jdata return np.array(self._eval) + +class PoyntingFlux(ObjectiveQuantity): + """A frequency-dependent Poynting Flux adjoint source. + Attributes: + volume: The volume over which the Poynting Flux is calculated. + This function currently only works for monitors with a defined + normal vector (e.g. planes in 3d or lines in 2d). User supplied + normal vectors may be implemented in the future. It also only + works with monitors aligned to a coordinate direction. + decimation_factor: Whether to skip points in the time series every + decimation_factor timesteps. See "add_dft_fields" documentation. + The biggest warning there is to be careful to avoid aliasing if + the fields vary quickly in time. + Note on yee_grid: For the Poynting Flux to work, H and E components + must lie at the same points. Therefore, the Yee grid will always be false. + """ + + def __init__(self, sim,volume, decimation_factor=0): + super().__init__(sim) + self.volume = sim._fit_volume_to_simulation(volume) + self.decimation_factor = decimation_factor + #_fit_volume_to_simulation increases the dimensionality of + #the volume, so we'll use the user input volume + self.normal = self.get_normal(volume) + + #Slightly different from other objectives, this + #function returns a list of monitors + def register_monitors(self, frequencies): + #Some template vectors to take the cross product with + #unfortunately Ex is 0 which I want to mean "no component", + #so I'm adding a fraction then rounding it out when passing + #which component fourierfields will grab + E_components = mp.Vector3(mp.Ex+0.1,mp.Ey,mp.Ez) + H_components = mp.Vector3(mp.Hx,mp.Hy,mp.Hz) + + #These are the 4 components of interest transverse + #to the plane, two for electric, two for magnetic + J_source_comps = H_components.cross(self.normal) + K_source_comps = self.normal.cross(E_components) + self.F_fields_list = [] + self.monitor_list = [] + #create FourierFields monitors for the electric currents + #(for the Poynting Flux adjoint, Ja = 1/4 H x n) + #Js is the component that FourierFields will retrieve + #e.g. two of Hx, Hy, or Hz + for Js in J_source_comps: + #There should always be at least one zero + if Js != 0: + self.F_fields_list + F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Js)), scale = (-1/4 if Js <0 else 1/4)) + self.F_fields_list.append(F_fields_here) + self.monitor_list.append(F_fields_here.register_monitors(frequencies)) + + #create FourierFields monitors for the magnetic currents + #(for the Poynting Flux adjoint, Ka = 1/4 n x E) + #Ks is the component that FourierFields will retrieve + #e.g. two of Ex, Ey, or Ez + for Ks in K_source_comps: + #There should always be at least one zero + if Ks != 0: + self.F_fields_list + F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Ks)), scale = (-1/4 if Ks <0 else 1/4)) + self.F_fields_list.append(F_fields_here) + self.monitor_list.append(F_fields_here.register_monitors(frequencies)) + + return self.monitor_list + + + def place_adjoint_source(self, dJ): + sources = [] + + for F_fields in self.F_fields_list: + sources.append(F_fields.place_adjoint_source(dJ)) + return sources + + def get_normal(self,volume): + #I'll add cylindrical later (since the normal vector gets a little different) + if (volume.dims == 2): + if (volume.size.x == 0): + return mp.Vector3(1,0,0) + elif(volume.size.y == 0): + return mp.Vector3(0,1,0) + else: + return mp.Vector3(0,0,1) + elif (volume.dims ==1): + if (volume.size.x == 0): + return mp.Vector3(1,0) + else: + return mp.Vector3(0,1) + + else : + return "Supported volumes are 1d or 2d" + + #This returns an array containing the Poynting Flux at each frequency point + #Note the first two FourierFields objectives are H fields, the second + #two are E fields. + #This covers the three cases of normal vectors using if statements. + #This assumes all the normal vectors point in the positive x,y,or z + #direction. The user may need to multiply by -1 to get the direction + #they expect if they're doing, for example, a box. + def __call__(self): + #it turns out the normal vector is the same equation in the x and z directions, but needs a negative 1 in the y direction + self._eval =((np.real(self.F_fields_list[0]()*np.conj(self.F_fields_list[3]()) - self.F_fields_list[1]()*np.conj(self.F_fields_list[2]())))) + if(self.normal.y == 1): + self._eval = -self.eval + + return self._eval + diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 56cee371a..3a63c34ca 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -13,7 +13,7 @@ from autograd import tensor_jacobian_product from utils import ApproxComparisonTestCase -MonitorObject = Enum("MonitorObject", "EIGENMODE DFT LDOS") +MonitorObject = Enum("MonitorObject", "EIGENMODE DFT LDOS POYNTING") class TestAdjointSolver(ApproxComparisonTestCase): @@ -169,6 +169,17 @@ def J(mode_mon): def J(mode_mon): return npa.array(mode_mon) + + elif mon_type.name == "POYNTING": + obj_list = [ + mpa.PoyntingFlux( + sim, + mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0.25, 1, 0)) + ) + ] + + def J(mode_mon): + return npa.sum(mode_mon) opt = mpa.OptimizationProblem( simulation=sim, @@ -358,6 +369,34 @@ def test_DFT_fields(self): tol = 0.07 if mp.is_single_precision() else 0.006 self.assertClose(adj_dd, fnd_dd, epsilon=tol) + + def test_Poynting_Flux(self): + print("*** TESTING POYNTING OBJECTIVE ***") + + # test the single frequency and multi frequency cases + for frequencies in [[self.fcen], [1 / 1.58, self.fcen, 1 / 1.53]]: + # compute objective value and its gradient for unperturbed design + unperturbed_val, unperturbed_grad = self.adjoint_solver( + self.p, MonitorObject.POYNTING, frequencies + ) + + # compute objective value for perturbed design + perturbed_val = self.adjoint_solver( + self.p + self.dp, MonitorObject.POYNTING, frequencies, need_gradient=False + ) + + # compare directional derivative + if unperturbed_grad.ndim < 2: + unperturbed_grad = np.expand_dims(unperturbed_grad, axis=1) + adj_dd = (self.dp[None, :] @ unperturbed_grad).flatten() + fnd_dd = perturbed_val - unperturbed_val + print( + f"directional derivative:, {adj_dd} (adjoint solver), {fnd_dd} (finite difference)" + ) + + tol = 0.07 if mp.is_single_precision() else 0.006 + self.assertClose(adj_dd, fnd_dd, epsilon=tol) + def test_eigenmode(self): print("*** TESTING EIGENMODE OBJECTIVE ***") From a0b1e0e5f7cbeb1b8b239b481ff7e11ce26585bb Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Fri, 5 Aug 2022 12:32:05 -0500 Subject: [PATCH 02/10] Added test file test_poynting_flux.py in main directory to reproduce bugs. --- test_poynting_flux.py | 132 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 test_poynting_flux.py diff --git a/test_poynting_flux.py b/test_poynting_flux.py new file mode 100644 index 000000000..9f44ef324 --- /dev/null +++ b/test_poynting_flux.py @@ -0,0 +1,132 @@ +import meep as mp +import numpy as np +import meep.adjoint as mpa +from autograd import numpy as npa +from autograd import tensor_jacobian_product +import nlopt +from matplotlib import pyplot as plt + +resolution = 30 # pixels/μm + +silicon = mp.Medium(epsilon=12) +sapphire = mp.Medium(epsilon_diag=(10.225,10.225,9.95), + epsilon_offdiag=(-0.825,-0.55*np.sqrt(3/2),0.55*np.sqrt(3/2))) + +sxy = 5.0 +cell_size = mp.Vector3(sxy,sxy,0) +#mp.integrate_field_function() + +dpml = 1.0 +pml_xy = [mp.PML(thickness=dpml)] +pml_x = [mp.PML(thickness=dpml,direction=mp.X)] + +eig_parity = mp.EVEN_Y+mp.ODD_Z + +design_region_size = mp.Vector3(1.5,1.5) +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(9861548) + +# random design region +p = 0.5*rng.rand(Nx*Ny) + +# random perturbation for design region +deps = 1e-5 +dp = deps*rng.rand(Nx*Ny) + +w = 1.0 +waveguide_geometry = [mp.Block(material=silicon, + center=mp.Vector3(), + size=mp.Vector3(mp.inf,w,mp.inf))] + +fcen = 1/1.55 +df = 0.2*fcen +mode_source = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df), + center=mp.Vector3(-0.5*sxy+dpml,0), + size=mp.Vector3(0,sxy-2*dpml), + eig_parity=eig_parity)] + +pt_source = [mp.Source(src=mp.GaussianSource(fcen,fwidth=df), + center=mp.Vector3(-0.5*sxy+dpml,0), + size=mp.Vector3(), + component=mp.Ez)] + +line_source = [mp.Source(src=mp.GaussianSource(fcen,fwidth=df), + center=mp.Vector3(-0.85,0), + size=mp.Vector3(0,sxy-2*dpml), + component=mp.Ez)] + +k_point = mp.Vector3(0.23,-0.38) + + +matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny), + mp.air, + silicon, + weights=np.ones((Nx,Ny))) + +matgrid_region = mpa.DesignRegion(matgrid, + volume=mp.Volume(center=mp.Vector3(), + size=mp.Vector3(design_region_size.x, + design_region_size.y, + 0))) + +matgrid_geometry = [mp.Block(center=matgrid_region.center, + size=matgrid_region.size, + material=matgrid)] + +geometry = waveguide_geometry + matgrid_geometry + +sim = mp.Simulation(resolution=resolution, + cell_size=cell_size, + boundary_layers=pml_xy, + sources=mode_source, + geometry=geometry) + + +frequencies = [fcen] +volume = mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0,1,0)) +obj_list = [mpa.PoyntingFlux(sim, volume)] +eval_hist = [] +def J(mode_mon): + eval = npa.sum(mode_mon) + return eval +opt = mpa.OptimizationProblem(simulation=sim, + objective_functions=J, + objective_arguments=obj_list, + design_regions=[matgrid_region], + frequencies=frequencies) + +f, dJ_du = opt([p]) +algorithm = nlopt.LD_MMA +n = Nx*Ny +maxeval = 10 + +# #all_fouriersrcdata = monitor.swigobj.fourier_sourcedata(self.volume.swigobj, self.component, self.sim.fields, dJ) + +# evaluation_history = [] +# sensitivity = [0] +# def f_fun(x, grad): +# f0, dJ_du = opt([x]) +# f0 = f0[0] # f0 is an array of length 1 +# if grad.size > 0: +# grad[:] = np.squeeze(dJ_du) +# evaluation_history.append(np.real(f0)) +# sensitivity[0] = dJ_du +# return np.real(f0) + + +# solver = nlopt.opt(algorithm, n) +# solver.set_lower_bounds(0) +# solver.set_upper_bounds(1) +# solver.set_max_objective(f_fun) +# solver.set_maxeval(maxeval) +# x = solver.optimize(p) +# plt.figure() +# plt.plot(evaluation_history,'o-') +# plt.grid(True) +# plt.xlabel('Iteration') +# plt.ylabel('FOM') +# plt.show() \ No newline at end of file From 2bf3eb51e61f09b5b66f129f934e02b82da2cc1c Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Sat, 6 Aug 2022 14:24:55 -0500 Subject: [PATCH 03/10] Fixed forward evalution to match the return of meep's poynting flux function. The Poynting Flux adjoint still fails when given field components of zero. The test file test_poynting_flux.py does the same simulation for adjoint and non-adjoint calls now. I'll need to add it into the actual test cases next, however. --- python/adjoint/objective.py | 33 +++++++++++++++++++++------- test_poynting_flux.py | 43 +++++++++++++++++++++++++++++++------ 2 files changed, 62 insertions(+), 14 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 1ae128327..3cc51785d 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -10,8 +10,7 @@ from .filter_source import FilteredSource Grid = namedtuple("Grid", ["x", "y", "z", "w"]) -E_CPTS = [mp.Ex, mp.Ey, mp.Ez] -H_CPTS = [mp.Hx, mp.Hy, mp.Hz] + class ObjectiveQuantity(abc.ABC): @@ -324,7 +323,12 @@ def place_adjoint_source(self, dJ): for fourier_data in self.all_fouriersrcdata: amp_arr = np.array(fourier_data.amp_arr).reshape(-1, self.num_freq) scale = self.scale * amp_arr * self._adj_src_scale(include_resolution=False) - + print("HERE'S THE SCALE:") + print("---------------------------") + print(scale) + print("----------------------------") + print() + print() if self.conj: scale = np.conjugate(scale) @@ -539,7 +543,7 @@ def register_monitors(self, frequencies): #There should always be at least one zero if Js != 0: self.F_fields_list - F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Js)), scale = (-1/4 if Js <0 else 1/4)) + F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Js)), scale = 1)#(-1/4 if Js <0 else 1/4)) self.F_fields_list.append(F_fields_here) self.monitor_list.append(F_fields_here.register_monitors(frequencies)) @@ -551,7 +555,7 @@ def register_monitors(self, frequencies): #There should always be at least one zero if Ks != 0: self.F_fields_list - F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Ks)), scale = (-1/4 if Ks <0 else 1/4)) + F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Ks)), scale =1)# (-1/4 if Ks <0 else 1/4)) self.F_fields_list.append(F_fields_here) self.monitor_list.append(F_fields_here.register_monitors(frequencies)) @@ -560,9 +564,10 @@ def register_monitors(self, frequencies): def place_adjoint_source(self, dJ): sources = [] - + i = 0 for F_fields in self.F_fields_list: - sources.append(F_fields.place_adjoint_source(dJ)) + if (np.abs(self.field_component_evaluations[i]) != 0): + sources.append(F_fields.place_adjoint_source(dJ)) return sources def get_normal(self,volume): @@ -592,7 +597,19 @@ def get_normal(self,volume): #they expect if they're doing, for example, a box. def __call__(self): #it turns out the normal vector is the same equation in the x and z directions, but needs a negative 1 in the y direction - self._eval =((np.real(self.F_fields_list[0]()*np.conj(self.F_fields_list[3]()) - self.F_fields_list[1]()*np.conj(self.F_fields_list[2]())))) + self.field_component_evaluations = [] + for field in self.F_fields_list: + field_here = field() + print(field_here) + if (np.size(field_here)==1): + field_here = np.zeros(1)#np.shape(1)) + self.field_component_evaluations.append(field_here) + integral_volume = (1 if self.volume.size.x == 0 else self.volume.size.x)*(1 if self.volume.size.y == 0 else self.volume.size.y)*(1 if self.volume.size.z == 0 else self.volume.size.z) + + #subtracting one makes it closer to meep's flux (probably some sort of array size thing) + discretization_factor = (self.field_component_evaluations[1].size/integral_volume)-1 + print(discretization_factor) + self._eval =np.sum((self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3]) - np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1]))/discretization_factor if(self.normal.y == 1): self._eval = -self.eval diff --git a/test_poynting_flux.py b/test_poynting_flux.py index 9f44ef324..b0d9ef88a 100644 --- a/test_poynting_flux.py +++ b/test_poynting_flux.py @@ -31,7 +31,8 @@ rng = np.random.RandomState(9861548) # random design region -p = 0.5*rng.rand(Nx*Ny) +#p = 0.5*rng.rand(Nx*Ny) +p = np.ones(Nx*Ny) # random perturbation for design region deps = 1e-5 @@ -87,7 +88,35 @@ frequencies = [fcen] -volume = mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0,1,0)) +print(fcen) +volume = mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0,2,0)) +# interesting_components = [mp.Hz,mp.Hy,mp.Ez,mp.Ey] +# dft_fields = sim.add_dft_fields(interesting_components, frequencies, where = volume,yee_grid = False) +# sim.run(until_after_sources=mp.stop_when_fields_decayed(100,mp.Ez,mp.Vector3(1.25,0),1e-12)) +# dft_field_data = [] +# for comp in interesting_components: +# dft_field_data.append(sim.get_dft_array(dft_fields,comp,0)) +# print(dft_field_data) +# discretization_factor = 1/2*dft_field_data[1].size +# print(discretization_factor) +# manual_flux_from_dft_fields = -np.sum(np.conj(dft_field_data[2])*dft_field_data[1])/(discretization_factor) +# print("manual_flux_from_dft_fields:") +# print(manual_flux_from_dft_fields) +# manual_flux = np.conj(dft_field_data) + +# refl_fr = mp.FluxRegion(center=mp.Vector3(1.25,0), size=mp.Vector3(0,2,0)) +# refl = sim.add_flux(fcen, df, 1,refl_fr) +# sim.run(until_after_sources=mp.stop_when_fields_decayed(100,mp.Ez,mp.Vector3(1.25,0),1e-2)) + +# straight_tran_data = sim.get_flux_data(refl) +# manual_flux = np.conj(straight_tran_data.H)*(straight_tran_data.E) +# print("manual flux is:") +# print(np.sum(manual_flux)) +# straight_tran_flux = mp.get_fluxes(refl) +# print(mp.get_flux_freqs(refl)) +# print(straight_tran_data) +# print(straight_tran_flux) + obj_list = [mpa.PoyntingFlux(sim, volume)] eval_hist = [] def J(mode_mon): @@ -99,10 +128,12 @@ def J(mode_mon): design_regions=[matgrid_region], frequencies=frequencies) -f, dJ_du = opt([p]) -algorithm = nlopt.LD_MMA -n = Nx*Ny -maxeval = 10 +f, dJ_du = opt([p],need_gradient=False) +print(f) + +# algorithm = nlopt.LD_MMA +# n = Nx*Ny +# maxeval = 10 # #all_fouriersrcdata = monitor.swigobj.fourier_sourcedata(self.volume.swigobj, self.component, self.sim.fields, dJ) From e3976e9d87735a3693d13c2a560714b7f7cdc2ce Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Tue, 9 Aug 2022 13:40:53 -0500 Subject: [PATCH 04/10] Forward run works using dft_fields objects. Adjoint run throws an error in step_db.cpp - Assertion 'changed_materials' failed. Creating a pull request will let me difftool code changes to debug the problem. --- python/adjoint/objective.py | 167 ++++++++++++++-------------- python/tests/test_adjoint_solver.py | 4 +- test_poynting_flux.py | 163 --------------------------- 3 files changed, 87 insertions(+), 247 deletions(-) delete mode 100644 test_poynting_flux.py diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 3cc51785d..fa5104214 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -11,7 +11,17 @@ Grid = namedtuple("Grid", ["x", "y", "z", "w"]) - +#3 possible components for E x n and H x n +#signs are handled in code +EH_TRANSVERSE = [ [mp.Hy, mp.Hz, mp.Ey, mp.Ez], + [mp.Hz, mp.Hx, mp.Ez, mp.Ex], + [mp.Hx, mp.Hy, mp.Ex, mp.Ey] ] + +#ToDO fix other two +#and comment too +JK_TRANSVERSE = [ [mp.Ey, mp.Hz, mp.Ex, mp.Hy], + [mp.Hz, mp.Hx, mp.Ez, mp.Ex], + [mp.Hx, mp.Hy, mp.Ex, mp.Ey] ] class ObjectiveQuantity(abc.ABC): """A differentiable objective quantity. @@ -269,14 +279,12 @@ def __call__(self): class FourierFields(ObjectiveQuantity): - def __init__(self, sim, volume, component, yee_grid=False, decimation_factor=0,scale = 1,conj = False): + def __init__(self, sim, volume, component, yee_grid=False, decimation_factor=0): super().__init__(sim) self.volume = sim._fit_volume_to_simulation(volume) self.component = component self.yee_grid = yee_grid self.decimation_factor = decimation_factor - self.scale = scale - self.conj = conj def register_monitors(self, frequencies): self._frequencies = np.asarray(frequencies) @@ -322,17 +330,10 @@ def place_adjoint_source(self, dJ): for fourier_data in self.all_fouriersrcdata: amp_arr = np.array(fourier_data.amp_arr).reshape(-1, self.num_freq) - scale = self.scale * amp_arr * self._adj_src_scale(include_resolution=False) - print("HERE'S THE SCALE:") - print("---------------------------") - print(scale) - print("----------------------------") - print() - print() - if self.conj: - scale = np.conjugate(scale) + scale = amp_arr * self._adj_src_scale(include_resolution=False) if self.num_freq == 1: + print(fourier_data.this) sources += [ mp.IndexedSource( time_src, fourier_data, scale[:, 0], not self.yee_grid @@ -513,105 +514,107 @@ class PoyntingFlux(ObjectiveQuantity): def __init__(self, sim,volume, decimation_factor=0): super().__init__(sim) - self.volume = sim._fit_volume_to_simulation(volume) - self.decimation_factor = decimation_factor #_fit_volume_to_simulation increases the dimensionality of #the volume, so we'll use the user input volume + self.volume = sim._fit_volume_to_simulation(volume) + self.decimation_factor = decimation_factor + #get_normal returns an index for the two + # dictionaries of cross products self.normal = self.get_normal(volume) + - #Slightly different from other objectives, this - #function returns a list of monitors + def register_monitors(self, frequencies): - #Some template vectors to take the cross product with - #unfortunately Ex is 0 which I want to mean "no component", - #so I'm adding a fraction then rounding it out when passing - #which component fourierfields will grab - E_components = mp.Vector3(mp.Ex+0.1,mp.Ey,mp.Ez) - H_components = mp.Vector3(mp.Hx,mp.Hy,mp.Hz) - - #These are the 4 components of interest transverse - #to the plane, two for electric, two for magnetic - J_source_comps = H_components.cross(self.normal) - K_source_comps = self.normal.cross(E_components) - self.F_fields_list = [] - self.monitor_list = [] - #create FourierFields monitors for the electric currents - #(for the Poynting Flux adjoint, Ja = 1/4 H x n) - #Js is the component that FourierFields will retrieve - #e.g. two of Hx, Hy, or Hz - for Js in J_source_comps: - #There should always be at least one zero - if Js != 0: - self.F_fields_list - F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Js)), scale = 1)#(-1/4 if Js <0 else 1/4)) - self.F_fields_list.append(F_fields_here) - self.monitor_list.append(F_fields_here.register_monitors(frequencies)) + self._frequencies = np.asarray(frequencies) - #create FourierFields monitors for the magnetic currents - #(for the Poynting Flux adjoint, Ka = 1/4 n x E) - #Ks is the component that FourierFields will retrieve - #e.g. two of Ex, Ey, or Ez - for Ks in K_source_comps: - #There should always be at least one zero - if Ks != 0: - self.F_fields_list - F_fields_here = FourierFields(self.sim,self.volume,np.around(abs(Ks)), scale =1)# (-1/4 if Ks <0 else 1/4)) - self.F_fields_list.append(F_fields_here) - self.monitor_list.append(F_fields_here.register_monitors(frequencies)) + #Add the dft monitors for the interesting components + self._monitor =self.sim.add_dft_fields(EH_TRANSVERSE[self.normal], + frequencies, + where = self.volume, + yee_grid = False) - return self.monitor_list + return self._monitor def place_adjoint_source(self, dJ): - sources = [] - i = 0 - for F_fields in self.F_fields_list: - if (np.abs(self.field_component_evaluations[i]) != 0): - sources.append(F_fields.place_adjoint_source(dJ)) - return sources + dJ = np.atleast_1d(dJ) + if dJ.ndim == 2: + dJ = np.sum(dJ, axis=1) + time_src = self._create_time_profile() + scale = self._adj_src_scale() + if self._frequencies.size == 1: + amp = dJ * scale + src = time_src + else: + scale = dJ * scale + src = FilteredSource( + time_src.frequency, + self._frequencies, + scale, + self.sim.fields.dt, + ) + amp = 1 + source = [] + #TODO add self.norm to be the right index 0,1, or 2 + for field_t,field in enumerate(JK_TRANSVERSE[self.normal]): + #dft_fields returns a scalar value for components that aren't evaluated + #in these case, we don't need to add an adjoint source + if(self.field_component_evaluations[field_t].size != 1): + tuple_elements = np.zeros((1,self.metadata[3].size,1), dtype=np.complex128) + #TODO: Load up the source data for other normal vectors + #TODO: Remove this unnecessary for loop + for index,element in enumerate(self.field_component_evaluations[field_t]): + tuple_elements[0,index,0] = element + source.append( mp.Source( + src, + component = mp.Ex, + amplitude=amp, + size=self.volume.size, + center=self.volume.center, + amp_data = tuple_elements + )) + + return source + #returns 0,1, or 2 corresponding to x,y, or z normal vectors + #TODO: Handle user-specified normal vectors and cases when 2d + #is projected into a dimension other than z def get_normal(self,volume): #I'll add cylindrical later (since the normal vector gets a little different) if (volume.dims == 2): if (volume.size.x == 0): - return mp.Vector3(1,0,0) + return 0 elif(volume.size.y == 0): - return mp.Vector3(0,1,0) + return 1 else: - return mp.Vector3(0,0,1) + return 2 elif (volume.dims ==1): if (volume.size.x == 0): - return mp.Vector3(1,0) + return 0 else: - return mp.Vector3(0,1) + return 1 else : return "Supported volumes are 1d or 2d" - #This returns an array containing the Poynting Flux at each frequency point + #Returns an array containing the Poynting Flux at each frequency point #Note the first two FourierFields objectives are H fields, the second #two are E fields. - #This covers the three cases of normal vectors using if statements. #This assumes all the normal vectors point in the positive x,y,or z #direction. The user may need to multiply by -1 to get the direction - #they expect if they're doing, for example, a box. + #they expect if they're doing, for example, a box. + #TODO: Add an extra parameter to let the user negate the flux, similar + #to the default meep ponynting flux def __call__(self): - #it turns out the normal vector is the same equation in the x and z directions, but needs a negative 1 in the y direction self.field_component_evaluations = [] - for field in self.F_fields_list: - field_here = field() - print(field_here) - if (np.size(field_here)==1): - field_here = np.zeros(1)#np.shape(1)) + #Loop over the relevant field components for this case of normal vector + for field in EH_TRANSVERSE[self.normal]: + field_here = self.sim.get_dft_array(self.monitor_list,field,0) self.field_component_evaluations.append(field_here) - integral_volume = (1 if self.volume.size.x == 0 else self.volume.size.x)*(1 if self.volume.size.y == 0 else self.volume.size.y)*(1 if self.volume.size.z == 0 else self.volume.size.z) - - #subtracting one makes it closer to meep's flux (probably some sort of array size thing) - discretization_factor = (self.field_component_evaluations[1].size/integral_volume)-1 - print(discretization_factor) - self._eval =np.sum((self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3]) - np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1]))/discretization_factor - if(self.normal.y == 1): - self._eval = -self.eval - + #Get weights for integration from cubature rule + self.metadata = self.sim.get_array_metadata(dft_cell = self.monitor_list) + flux_density =np.real(-(self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3])) + (np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1])) + #Apply cubature weights + self._eval = np.sum(self.metadata[3]*flux_density) return self._eval diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 3a63c34ca..25dbca662 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -174,12 +174,12 @@ def J(mode_mon): obj_list = [ mpa.PoyntingFlux( sim, - mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0.25, 1, 0)) + mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0, 1, 0)) ) ] def J(mode_mon): - return npa.sum(mode_mon) + return mode_mon opt = mpa.OptimizationProblem( simulation=sim, diff --git a/test_poynting_flux.py b/test_poynting_flux.py deleted file mode 100644 index b0d9ef88a..000000000 --- a/test_poynting_flux.py +++ /dev/null @@ -1,163 +0,0 @@ -import meep as mp -import numpy as np -import meep.adjoint as mpa -from autograd import numpy as npa -from autograd import tensor_jacobian_product -import nlopt -from matplotlib import pyplot as plt - -resolution = 30 # pixels/μm - -silicon = mp.Medium(epsilon=12) -sapphire = mp.Medium(epsilon_diag=(10.225,10.225,9.95), - epsilon_offdiag=(-0.825,-0.55*np.sqrt(3/2),0.55*np.sqrt(3/2))) - -sxy = 5.0 -cell_size = mp.Vector3(sxy,sxy,0) -#mp.integrate_field_function() - -dpml = 1.0 -pml_xy = [mp.PML(thickness=dpml)] -pml_x = [mp.PML(thickness=dpml,direction=mp.X)] - -eig_parity = mp.EVEN_Y+mp.ODD_Z - -design_region_size = mp.Vector3(1.5,1.5) -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(9861548) - -# random design region -#p = 0.5*rng.rand(Nx*Ny) -p = np.ones(Nx*Ny) - -# random perturbation for design region -deps = 1e-5 -dp = deps*rng.rand(Nx*Ny) - -w = 1.0 -waveguide_geometry = [mp.Block(material=silicon, - center=mp.Vector3(), - size=mp.Vector3(mp.inf,w,mp.inf))] - -fcen = 1/1.55 -df = 0.2*fcen -mode_source = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df), - center=mp.Vector3(-0.5*sxy+dpml,0), - size=mp.Vector3(0,sxy-2*dpml), - eig_parity=eig_parity)] - -pt_source = [mp.Source(src=mp.GaussianSource(fcen,fwidth=df), - center=mp.Vector3(-0.5*sxy+dpml,0), - size=mp.Vector3(), - component=mp.Ez)] - -line_source = [mp.Source(src=mp.GaussianSource(fcen,fwidth=df), - center=mp.Vector3(-0.85,0), - size=mp.Vector3(0,sxy-2*dpml), - component=mp.Ez)] - -k_point = mp.Vector3(0.23,-0.38) - - -matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny), - mp.air, - silicon, - weights=np.ones((Nx,Ny))) - -matgrid_region = mpa.DesignRegion(matgrid, - volume=mp.Volume(center=mp.Vector3(), - size=mp.Vector3(design_region_size.x, - design_region_size.y, - 0))) - -matgrid_geometry = [mp.Block(center=matgrid_region.center, - size=matgrid_region.size, - material=matgrid)] - -geometry = waveguide_geometry + matgrid_geometry - -sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, - boundary_layers=pml_xy, - sources=mode_source, - geometry=geometry) - - -frequencies = [fcen] -print(fcen) -volume = mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0,2,0)) -# interesting_components = [mp.Hz,mp.Hy,mp.Ez,mp.Ey] -# dft_fields = sim.add_dft_fields(interesting_components, frequencies, where = volume,yee_grid = False) -# sim.run(until_after_sources=mp.stop_when_fields_decayed(100,mp.Ez,mp.Vector3(1.25,0),1e-12)) -# dft_field_data = [] -# for comp in interesting_components: -# dft_field_data.append(sim.get_dft_array(dft_fields,comp,0)) -# print(dft_field_data) -# discretization_factor = 1/2*dft_field_data[1].size -# print(discretization_factor) -# manual_flux_from_dft_fields = -np.sum(np.conj(dft_field_data[2])*dft_field_data[1])/(discretization_factor) -# print("manual_flux_from_dft_fields:") -# print(manual_flux_from_dft_fields) -# manual_flux = np.conj(dft_field_data) - -# refl_fr = mp.FluxRegion(center=mp.Vector3(1.25,0), size=mp.Vector3(0,2,0)) -# refl = sim.add_flux(fcen, df, 1,refl_fr) -# sim.run(until_after_sources=mp.stop_when_fields_decayed(100,mp.Ez,mp.Vector3(1.25,0),1e-2)) - -# straight_tran_data = sim.get_flux_data(refl) -# manual_flux = np.conj(straight_tran_data.H)*(straight_tran_data.E) -# print("manual flux is:") -# print(np.sum(manual_flux)) -# straight_tran_flux = mp.get_fluxes(refl) -# print(mp.get_flux_freqs(refl)) -# print(straight_tran_data) -# print(straight_tran_flux) - -obj_list = [mpa.PoyntingFlux(sim, volume)] -eval_hist = [] -def J(mode_mon): - eval = npa.sum(mode_mon) - return eval -opt = mpa.OptimizationProblem(simulation=sim, - objective_functions=J, - objective_arguments=obj_list, - design_regions=[matgrid_region], - frequencies=frequencies) - -f, dJ_du = opt([p],need_gradient=False) -print(f) - -# algorithm = nlopt.LD_MMA -# n = Nx*Ny -# maxeval = 10 - -# #all_fouriersrcdata = monitor.swigobj.fourier_sourcedata(self.volume.swigobj, self.component, self.sim.fields, dJ) - -# evaluation_history = [] -# sensitivity = [0] -# def f_fun(x, grad): -# f0, dJ_du = opt([x]) -# f0 = f0[0] # f0 is an array of length 1 -# if grad.size > 0: -# grad[:] = np.squeeze(dJ_du) -# evaluation_history.append(np.real(f0)) -# sensitivity[0] = dJ_du -# return np.real(f0) - - -# solver = nlopt.opt(algorithm, n) -# solver.set_lower_bounds(0) -# solver.set_upper_bounds(1) -# solver.set_max_objective(f_fun) -# solver.set_maxeval(maxeval) -# x = solver.optimize(p) -# plt.figure() -# plt.plot(evaluation_history,'o-') -# plt.grid(True) -# plt.xlabel('Iteration') -# plt.ylabel('FOM') -# plt.show() \ No newline at end of file From 425c07ab421d6d1d5ab974025c5b2a511388966a Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Tue, 9 Aug 2022 14:12:09 -0500 Subject: [PATCH 05/10] =Fixed JK_TRANSVERSE, removed some last edits in FourierFields, and added some more TODOs --- python/adjoint/objective.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index fa5104214..83af78158 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -17,11 +17,12 @@ [mp.Hz, mp.Hx, mp.Ez, mp.Ex], [mp.Hx, mp.Hy, mp.Ex, mp.Ey] ] -#ToDO fix other two -#and comment too -JK_TRANSVERSE = [ [mp.Ey, mp.Hz, mp.Ex, mp.Hy], - [mp.Hz, mp.Hx, mp.Ez, mp.Ex], - [mp.Hx, mp.Hy, mp.Ex, mp.Ey] ] +#Holds the components for each current source +#for the cases of x,y, and z normal vectors. +#This is the same as swapping H and E in the above list +JK_TRANSVERSE = [ [mp.Ey, mp.Ez, mp.Hy, mp.Hz], + [mp.Ez, mp.Ex, mp.Hz, mp.Hx], + [mp.Ex, mp.Ey, mp.Hx, mp.Hy] ] class ObjectiveQuantity(abc.ABC): """A differentiable objective quantity. @@ -333,7 +334,6 @@ def place_adjoint_source(self, dJ): scale = amp_arr * self._adj_src_scale(include_resolution=False) if self.num_freq == 1: - print(fourier_data.this) sources += [ mp.IndexedSource( time_src, fourier_data, scale[:, 0], not self.yee_grid @@ -555,10 +555,11 @@ def place_adjoint_source(self, dJ): ) amp = 1 source = [] - #TODO add self.norm to be the right index 0,1, or 2 + #TODO: account for sign in the current sources + #(the K sources should have the negative) for field_t,field in enumerate(JK_TRANSVERSE[self.normal]): #dft_fields returns a scalar value for components that aren't evaluated - #in these case, we don't need to add an adjoint source + #in these cases, we don't need to add an adjoint source if(self.field_component_evaluations[field_t].size != 1): tuple_elements = np.zeros((1,self.metadata[3].size,1), dtype=np.complex128) #TODO: Load up the source data for other normal vectors @@ -567,7 +568,7 @@ def place_adjoint_source(self, dJ): tuple_elements[0,index,0] = element source.append( mp.Source( src, - component = mp.Ex, + component = field, amplitude=amp, size=self.volume.size, center=self.volume.center, From bc121b78db2b061c21a5c9923332da8c347c8a9e Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Wed, 10 Aug 2022 13:00:03 -0500 Subject: [PATCH 06/10] Fixed JK_TRANSVERSE so the adjiont simulation finishes without error. Off by a factor of ten in initial tests, debugging today. --- python/adjoint/objective.py | 62 ++++++++++++++--------------- python/tests/test_adjoint_solver.py | 21 +++++----- 2 files changed, 39 insertions(+), 44 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 83af78158..0f551c9f5 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -20,9 +20,9 @@ #Holds the components for each current source #for the cases of x,y, and z normal vectors. #This is the same as swapping H and E in the above list -JK_TRANSVERSE = [ [mp.Ey, mp.Ez, mp.Hy, mp.Hz], - [mp.Ez, mp.Ex, mp.Hz, mp.Hx], - [mp.Ex, mp.Ey, mp.Hx, mp.Hy] ] +JK_TRANSVERSE = [ [mp.Ez, mp.Ey, mp.Hz, mp.Hy], + [mp.Ex, mp.Ez, mp.Hx, mp.Hz], + [mp.Ey, mp.Ex, mp.Hy, mp.Hx] ] class ObjectiveQuantity(abc.ABC): """A differentiable objective quantity. @@ -511,31 +511,34 @@ class PoyntingFlux(ObjectiveQuantity): Note on yee_grid: For the Poynting Flux to work, H and E components must lie at the same points. Therefore, the Yee grid will always be false. """ - - def __init__(self, sim,volume, decimation_factor=0): + def __init__( + self, + sim, + volume, + scale = 1, + decimation_factor=0 + ): super().__init__(sim) #_fit_volume_to_simulation increases the dimensionality of #the volume, so we'll use the user input volume self.volume = sim._fit_volume_to_simulation(volume) self.decimation_factor = decimation_factor + self.scale = scale #get_normal returns an index for the two # dictionaries of cross products self.normal = self.get_normal(volume) - - def register_monitors(self, frequencies): self._frequencies = np.asarray(frequencies) - + #Add the dft monitors for the interesting components self._monitor =self.sim.add_dft_fields(EH_TRANSVERSE[self.normal], frequencies, where = self.volume, yee_grid = False) - + return self._monitor - - + def place_adjoint_source(self, dJ): dJ = np.atleast_1d(dJ) if dJ.ndim == 2: @@ -574,8 +577,22 @@ def place_adjoint_source(self, dJ): center=self.volume.center, amp_data = tuple_elements )) - return source + + def __call__(self): + self.field_component_evaluations = [] + #Loop over the relevant field components for this case of normal vector + for field in EH_TRANSVERSE[self.normal]: + field_here = self.sim.get_dft_array(self._monitor,field,0) + self.field_component_evaluations.append(field_here) + #Get weights for integration from cubature rule + self.metadata = self.sim.get_array_metadata(dft_cell = self._monitor) + flux_density =np.real(-(self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3])) + (np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1])) + #Apply cubature weights and user-supplied scale + self._eval =self.scale * np.sum(self.metadata[3]*flux_density) + return self._eval + + #returns 0,1, or 2 corresponding to x,y, or z normal vectors #TODO: Handle user-specified normal vectors and cases when 2d @@ -598,24 +615,3 @@ def get_normal(self,volume): else : return "Supported volumes are 1d or 2d" - #Returns an array containing the Poynting Flux at each frequency point - #Note the first two FourierFields objectives are H fields, the second - #two are E fields. - #This assumes all the normal vectors point in the positive x,y,or z - #direction. The user may need to multiply by -1 to get the direction - #they expect if they're doing, for example, a box. - #TODO: Add an extra parameter to let the user negate the flux, similar - #to the default meep ponynting flux - def __call__(self): - self.field_component_evaluations = [] - #Loop over the relevant field components for this case of normal vector - for field in EH_TRANSVERSE[self.normal]: - field_here = self.sim.get_dft_array(self.monitor_list,field,0) - self.field_component_evaluations.append(field_here) - #Get weights for integration from cubature rule - self.metadata = self.sim.get_array_metadata(dft_cell = self.monitor_list) - flux_density =np.real(-(self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3])) + (np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1])) - #Apply cubature weights - self._eval = np.sum(self.metadata[3]*flux_density) - return self._eval - diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 25dbca662..485d3794b 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -342,20 +342,20 @@ def mapping(self, x, filter_radius, eta, beta): projected_field = mpa.tanh_projection(filtered_field, beta, eta) return projected_field.flatten() - - def test_DFT_fields(self): - print("*** TESTING DFT OBJECTIVE ***") + + def test_Poynting_Flux(self): + print("*** TESTING POYNTING OBJECTIVE ***") # test the single frequency and multi frequency cases for frequencies in [[self.fcen], [1 / 1.58, self.fcen, 1 / 1.53]]: # compute objective value and its gradient for unperturbed design unperturbed_val, unperturbed_grad = self.adjoint_solver( - self.p, MonitorObject.DFT, frequencies + self.p, MonitorObject.POYNTING, frequencies ) # compute objective value for perturbed design perturbed_val = self.adjoint_solver( - self.p + self.dp, MonitorObject.DFT, frequencies, need_gradient=False + self.p + self.dp, MonitorObject.POYNTING, frequencies, need_gradient=False ) # compare directional derivative @@ -369,20 +369,20 @@ def test_DFT_fields(self): tol = 0.07 if mp.is_single_precision() else 0.006 self.assertClose(adj_dd, fnd_dd, epsilon=tol) - - def test_Poynting_Flux(self): - print("*** TESTING POYNTING OBJECTIVE ***") + + def test_DFT_fields(self): + print("*** TESTING DFT OBJECTIVE ***") # test the single frequency and multi frequency cases for frequencies in [[self.fcen], [1 / 1.58, self.fcen, 1 / 1.53]]: # compute objective value and its gradient for unperturbed design unperturbed_val, unperturbed_grad = self.adjoint_solver( - self.p, MonitorObject.POYNTING, frequencies + self.p, MonitorObject.DFT, frequencies ) # compute objective value for perturbed design perturbed_val = self.adjoint_solver( - self.p + self.dp, MonitorObject.POYNTING, frequencies, need_gradient=False + self.p + self.dp, MonitorObject.DFT, frequencies, need_gradient=False ) # compare directional derivative @@ -397,7 +397,6 @@ def test_Poynting_Flux(self): tol = 0.07 if mp.is_single_precision() else 0.006 self.assertClose(adj_dd, fnd_dd, epsilon=tol) - def test_eigenmode(self): print("*** TESTING EIGENMODE OBJECTIVE ***") From 64a434b9abb0338971844adcf20f26c441e690c5 Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Tue, 16 Aug 2022 14:22:51 -0400 Subject: [PATCH 07/10] Fixed sign and conjugate bugs. Changed the implementation from four line sources (one for each tangential H or E component) to a series of point sources to accomodate end-to-end interpolation - we needed middle to middle. Unfortunately, this did not seem to improve the error appreciably. --- python/adjoint/objective.py | 99 ++++++++++++++++------------- python/tests/test_adjoint_solver.py | 42 ++++++++++++ 2 files changed, 98 insertions(+), 43 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 0f551c9f5..98956a697 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -13,16 +13,19 @@ #3 possible components for E x n and H x n #signs are handled in code -EH_TRANSVERSE = [ [mp.Hy, mp.Hz, mp.Ey, mp.Ez], - [mp.Hz, mp.Hx, mp.Ez, mp.Ex], - [mp.Hx, mp.Hy, mp.Ex, mp.Ey] ] +EH_TRANSVERSE = [ [mp.Hz, mp.Hy, mp.Ez, mp.Ey], + [mp.Hx, mp.Hz, mp.Ex, mp.Ez], + [mp.Hy, mp.Hx, mp.Ey, mp.Ex] ] #Holds the components for each current source #for the cases of x,y, and z normal vectors. #This is the same as swapping H and E in the above list -JK_TRANSVERSE = [ [mp.Ez, mp.Ey, mp.Hz, mp.Hy], - [mp.Ex, mp.Ez, mp.Hx, mp.Hz], - [mp.Ey, mp.Ex, mp.Hy, mp.Hx] ] +JK_TRANSVERSE = [ [mp.Ey, mp.Ez, mp.Hy, mp.Hz], + [mp.Ez, mp.Ex, mp.Hz, mp.Hx], + [mp.Ex, mp.Ey, mp.Hx, mp.Hy] ] + +#Holds the amplitudes used in Poynting Flux adjoint sources +FLUX_AMPLITUDES = np.array([1/4,-1/4,-1/4,1/4], dtype = np.complex128) class ObjectiveQuantity(abc.ABC): """A differentiable objective quantity. @@ -544,7 +547,7 @@ def place_adjoint_source(self, dJ): if dJ.ndim == 2: dJ = np.sum(dJ, axis=1) time_src = self._create_time_profile() - scale = self._adj_src_scale() + scale = self._adj_src_scale(include_resolution = False) if self._frequencies.size == 1: amp = dJ * scale src = time_src @@ -557,27 +560,42 @@ def place_adjoint_source(self, dJ): self.sim.fields.dt, ) amp = 1 + #Sometimes a couple field component arrays contain a single zero (typically + # due to symmetries). Make a list of where field components with multiple elements are + nonzero_comps = [] + for field_index,field_comp in enumerate(self.field_component_evaluations): + if(field_comp.size != 1): + nonzero_comps.append(field_index) source = [] - #TODO: account for sign in the current sources - #(the K sources should have the negative) - for field_t,field in enumerate(JK_TRANSVERSE[self.normal]): - #dft_fields returns a scalar value for components that aren't evaluated - #in these cases, we don't need to add an adjoint source - if(self.field_component_evaluations[field_t].size != 1): - tuple_elements = np.zeros((1,self.metadata[3].size,1), dtype=np.complex128) - #TODO: Load up the source data for other normal vectors - #TODO: Remove this unnecessary for loop - for index,element in enumerate(self.field_component_evaluations[field_t]): - tuple_elements[0,index,0] = element - source.append( mp.Source( - src, - component = field, - amplitude=amp, - size=self.volume.size, - center=self.volume.center, - amp_data = tuple_elements - )) + #Get the right components for this normal vector + field_components_for_sources = JK_TRANSVERSE[self.normal] + for idx in range(self.metadata[0].size): + for idy in range(self.metadata[1].size): + for idz in range(self.metadata[2].size): + #Iterate over the relevant indices in the constants defined above + for field_idx in nonzero_comps: + #dft_fields returns a scalar value for components that aren't evaluated + #in these cases, we don't need to add an adjoint source + field_component_array = self.field_component_evaluations[field_idx] + source_component = field_components_for_sources[field_idx] + #scale the amplitude by +/- 1/4 from FLUX_AMPLITUDES, then apply the amplitude at this spatial point + #from the field component + #TODO: use the not-yet-implemented uncollapse dimensions function to make field commponent array 3 dimensional + #It currently only works for x-directed normal vectors in a 2-d simulation + amplitude = FLUX_AMPLITUDES[field_idx]*np.conj(field_component_array[idy])*amp + center = mp.Vector3(self.metadata[0][idx],self.metadata[1][idy],self.metadata[2][idz]) + source.append(mp.Source(src, + component = source_component, + amplitude= amplitude, + center=center + )) return source + + #Will eventually turn 1xn arrays of component amplitudes along a boundary + #into 1xnxn or nx1x1, etc. arrays depending on the normal vector + #so that they can be looped over in three dimensions correctly + #def uncollapse_dimensions(self,component_array, metadata): + def __call__(self): self.field_component_evaluations = [] @@ -587,31 +605,26 @@ def __call__(self): self.field_component_evaluations.append(field_here) #Get weights for integration from cubature rule self.metadata = self.sim.get_array_metadata(dft_cell = self._monitor) - flux_density =np.real(-(self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3])) + (np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1])) + flux_density =np.real((self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3])) - (np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1])) #Apply cubature weights and user-supplied scale + # print("Here are the metadata weights:") + # print(self.metadata[3]) self._eval =self.scale * np.sum(self.metadata[3]*flux_density) + print("Here is the evaluation:") + print(self._eval) return self._eval - #returns 0,1, or 2 corresponding to x,y, or z normal vectors + #returns 0,1, or 2 corresponding to x, y, or z normal vectors #TODO: Handle user-specified normal vectors and cases when 2d #is projected into a dimension other than z def get_normal(self,volume): #I'll add cylindrical later (since the normal vector gets a little different) - if (volume.dims == 2): - if (volume.size.x == 0): - return 0 - elif(volume.size.y == 0): - return 1 - else: - return 2 - elif (volume.dims ==1): - if (volume.size.x == 0): - return 0 - else: - return 1 - - else : - return "Supported volumes are 1d or 2d" + if (volume.size.x == 0): + return 0 + elif(volume.size.y == 0): + return 1 + else: + return 2 diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 485d3794b..5f1935e71 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -366,9 +366,47 @@ def test_Poynting_Flux(self): print( f"directional derivative:, {adj_dd} (adjoint solver), {fnd_dd} (finite difference)" ) + print("\n\n\n\n\n") + print("") + print("This is the unperturbed poynting gradient:") + print(unperturbed_grad) tol = 0.07 if mp.is_single_precision() else 0.006 self.assertClose(adj_dd, fnd_dd, epsilon=tol) + + def test_eigenmode(self): + print("*** TESTING EIGENMODE OBJECTIVE ***") + + # test the single frequency and multi frequency case + for frequencies in [[self.fcen], [1 / 1.58, self.fcen, 1 / 1.53]]: + # compute objective value and its gradient for unperturbed design + unperturbed_val, unperturbed_grad = self.adjoint_solver( + self.p, MonitorObject.EIGENMODE, frequencies + ) + + # compute objective for perturbed design + perturbed_val = self.adjoint_solver( + self.p + self.dp, + MonitorObject.EIGENMODE, + frequencies, + need_gradient=False, + ) + + # compare directional derivative + if unperturbed_grad.ndim < 2: + unperturbed_grad = np.expand_dims(unperturbed_grad, axis=1) + adj_dd = (self.dp[None, :] @ unperturbed_grad).flatten() + fnd_dd = perturbed_val - unperturbed_val + print( + f"directional derivative:, {adj_dd} (adjoint solver), {fnd_dd} (finite difference)" + ) + print("\n\n\n\n\n") + print("") + print("This is the unperturbed eigenmode gradient:") + print(unperturbed_grad) + + tol = 0.04 if mp.is_single_precision() else 0.01 + self.assertClose(adj_dd, fnd_dd, epsilon=tol) def test_DFT_fields(self): print("*** TESTING DFT OBJECTIVE ***") @@ -423,6 +461,10 @@ def test_eigenmode(self): print( f"directional derivative:, {adj_dd} (adjoint solver), {fnd_dd} (finite difference)" ) + print("\n\n\n\n\n") + print("") + print("This is the unperturbed eigenmode gradient:") + print(unperturbed_grad) tol = 0.04 if mp.is_single_precision() else 0.01 self.assertClose(adj_dd, fnd_dd, epsilon=tol) From 475594a0a1525f8b1bdf119233fe189cf01cf550 Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Thu, 25 Aug 2022 11:55:09 -0400 Subject: [PATCH 08/10] Changed implementation back to FourierFields, passing the gradient as dJ. Works well in 2D, fails in 3D. Multiple frequencies untested. --- python/adjoint/objective.py | 187 ++++++++++--------------- python/adjoint/optimization_problem.py | 2 + python/tests/test_adjoint_solver.py | 30 ++-- 3 files changed, 96 insertions(+), 123 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 98956a697..54f94adeb 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -1,7 +1,7 @@ """Handling of objective functions and objective quantities.""" import abc from collections import namedtuple - +from autograd import numpy as npa import numpy as np from meep.simulation import py_v3_to_vec @@ -11,21 +11,26 @@ Grid = namedtuple("Grid", ["x", "y", "z", "w"]) -#3 possible components for E x n and H x n -#signs are handled in code -EH_TRANSVERSE = [ [mp.Hz, mp.Hy, mp.Ez, mp.Ey], - [mp.Hx, mp.Hz, mp.Ex, mp.Ez], - [mp.Hy, mp.Hx, mp.Ey, mp.Ex] ] +# 3 possible components for E x n and H x n +# signs are handled in code +EH_TRANSVERSE = [ + [mp.Hz, mp.Hy, mp.Ez, mp.Ey], + [mp.Hx, mp.Hz, mp.Ex, mp.Ez], + [mp.Hy, mp.Hx, mp.Ey, mp.Ex], +] + +# Holds the components for each current source +# for the cases of x,y, and z normal vectors. +# This is the same as swapping H and E in the above list +JK_TRANSVERSE = [ + [mp.Ey, mp.Ez, mp.Hy, mp.Hz], + [mp.Ez, mp.Ex, mp.Hz, mp.Hx], + [mp.Ex, mp.Ey, mp.Hx, mp.Hy], +] -#Holds the components for each current source -#for the cases of x,y, and z normal vectors. -#This is the same as swapping H and E in the above list -JK_TRANSVERSE = [ [mp.Ey, mp.Ez, mp.Hy, mp.Hz], - [mp.Ez, mp.Ex, mp.Hz, mp.Hx], - [mp.Ex, mp.Ey, mp.Hx, mp.Hy] ] +# Holds the amplitudes used in Poynting Flux adjoint sources +FLUX_AMPLITUDES = np.array([1 / 4, -1 / 4, -1 / 4, 1 / 4], dtype=np.complex128) -#Holds the amplitudes used in Poynting Flux adjoint sources -FLUX_AMPLITUDES = np.array([1/4,-1/4,-1/4,1/4], dtype = np.complex128) class ObjectiveQuantity(abc.ABC): """A differentiable objective quantity. @@ -331,7 +336,7 @@ def place_adjoint_source(self, dJ): self.all_fouriersrcdata = self._monitor.swigobj.fourier_sourcedata( self.volume.swigobj, self.component, self.sim.fields, dJ ) - + for fourier_data in self.all_fouriersrcdata: amp_arr = np.array(fourier_data.amp_arr).reshape(-1, self.num_freq) scale = amp_arr * self._adj_src_scale(include_resolution=False) @@ -498,7 +503,8 @@ def __call__(self): self.ldos_scale = self.sim.ldos_scale self.ldos_Jdata = self.sim.ldos_Jdata return np.array(self._eval) - + + class PoyntingFlux(ObjectiveQuantity): """A frequency-dependent Poynting Flux adjoint source. Attributes: @@ -507,124 +513,85 @@ class PoyntingFlux(ObjectiveQuantity): normal vector (e.g. planes in 3d or lines in 2d). User supplied normal vectors may be implemented in the future. It also only works with monitors aligned to a coordinate direction. - decimation_factor: Whether to skip points in the time series every + decimation_factor: Whether to skip points in the time series every decimation_factor timesteps. See "add_dft_fields" documentation. The biggest warning there is to be careful to avoid aliasing if the fields vary quickly in time. Note on yee_grid: For the Poynting Flux to work, H and E components must lie at the same points. Therefore, the Yee grid will always be false. """ - def __init__( - self, - sim, - volume, - scale = 1, - decimation_factor=0 - ): + + def __init__(self, sim, volume, scale=1, decimation_factor=0): super().__init__(sim) - #_fit_volume_to_simulation increases the dimensionality of - #the volume, so we'll use the user input volume + # _fit_volume_to_simulation increases the dimensionality of + # the volume, so we'll use the user input volume self.volume = sim._fit_volume_to_simulation(volume) self.decimation_factor = decimation_factor self.scale = scale - #get_normal returns an index for the two + # get_normal returns an index for the two # dictionaries of cross products self.normal = self.get_normal(volume) def register_monitors(self, frequencies): self._frequencies = np.asarray(frequencies) - - #Add the dft monitors for the interesting components - self._monitor =self.sim.add_dft_fields(EH_TRANSVERSE[self.normal], - frequencies, - where = self.volume, - yee_grid = False) - + self._monitor = [] + # List to hold FourierFields objects + self.F_fields_list = [] + for comp in EH_TRANSVERSE[self.normal]: + # instantiate the FourierFields monitors + F_field = FourierFields(self.sim, self.volume, comp) + self.F_fields_list.append(F_field) + self._monitor.append(F_field.register_monitors(self._frequencies)) return self._monitor def place_adjoint_source(self, dJ): - dJ = np.atleast_1d(dJ) - if dJ.ndim == 2: - dJ = np.sum(dJ, axis=1) - time_src = self._create_time_profile() - scale = self._adj_src_scale(include_resolution = False) - if self._frequencies.size == 1: - amp = dJ * scale - src = time_src - else: - scale = dJ * scale - src = FilteredSource( - time_src.frequency, - self._frequencies, - scale, - self.sim.fields.dt, - ) - amp = 1 - #Sometimes a couple field component arrays contain a single zero (typically - # due to symmetries). Make a list of where field components with multiple elements are - nonzero_comps = [] - for field_index,field_comp in enumerate(self.field_component_evaluations): - if(field_comp.size != 1): - nonzero_comps.append(field_index) source = [] - #Get the right components for this normal vector - field_components_for_sources = JK_TRANSVERSE[self.normal] - for idx in range(self.metadata[0].size): - for idy in range(self.metadata[1].size): - for idz in range(self.metadata[2].size): - #Iterate over the relevant indices in the constants defined above - for field_idx in nonzero_comps: - #dft_fields returns a scalar value for components that aren't evaluated - #in these cases, we don't need to add an adjoint source - field_component_array = self.field_component_evaluations[field_idx] - source_component = field_components_for_sources[field_idx] - #scale the amplitude by +/- 1/4 from FLUX_AMPLITUDES, then apply the amplitude at this spatial point - #from the field component - #TODO: use the not-yet-implemented uncollapse dimensions function to make field commponent array 3 dimensional - #It currently only works for x-directed normal vectors in a 2-d simulation - amplitude = FLUX_AMPLITUDES[field_idx]*np.conj(field_component_array[idy])*amp - center = mp.Vector3(self.metadata[0][idx],self.metadata[1][idy],self.metadata[2][idz]) - source.append(mp.Source(src, - component = source_component, - amplitude= amplitude, - center=center - )) + for pos, field in enumerate(self.F_fields_list): + source.append(field.place_adjoint_source(np.flipud(dJ[pos].flatten()))[0]) return source - - #Will eventually turn 1xn arrays of component amplitudes along a boundary - #into 1xnxn or nx1x1, etc. arrays depending on the normal vector - #so that they can be looped over in three dimensions correctly - #def uncollapse_dimensions(self,component_array, metadata): - def __call__(self): self.field_component_evaluations = [] - #Loop over the relevant field components for this case of normal vector - for field in EH_TRANSVERSE[self.normal]: - field_here = self.sim.get_dft_array(self._monitor,field,0) - self.field_component_evaluations.append(field_here) - #Get weights for integration from cubature rule - self.metadata = self.sim.get_array_metadata(dft_cell = self._monitor) - flux_density =np.real((self.field_component_evaluations[0]*np.conj(self.field_component_evaluations[3])) - (np.conj(self.field_component_evaluations[2])*self.field_component_evaluations[1])) - #Apply cubature weights and user-supplied scale - # print("Here are the metadata weights:") - # print(self.metadata[3]) - self._eval =self.scale * np.sum(self.metadata[3]*flux_density) - print("Here is the evaluation:") - print(self._eval) - return self._eval - - - - #returns 0,1, or 2 corresponding to x, y, or z normal vectors - #TODO: Handle user-specified normal vectors and cases when 2d - #is projected into a dimension other than z - def get_normal(self,volume): - #I'll add cylindrical later (since the normal vector gets a little different) - if (volume.size.x == 0): + for field in self.F_fields_list: + # Get the dft evaluation from a call to the underlying + # FourierFields object + field_here = field() + self.field_component_evaluations.append(field_here) + # Get integration weights Meep uses + self.metadata = self.sim.get_array_metadata(vol=self.volume) + self.field_component_evaluations.append( + np.array([self.metadata[3]]).astype(complex) + ) + [H1, H2, E1, E2, meta] = self.field_component_evaluations + self._eval = self.field_component_evaluations + print("This is meta*E2") + print(meta * E2) + return np.array([H1, H2, E1, E2, meta]) + + # takes in a 1x5xNxM vector where the size five array corresponds to + # [H1,H2,E1,E1,meta] + # multiple frequencies will be tested later + @staticmethod + def compute_flux(*inputs): + flux = npa.sum( + npa.real( + inputs[0][4] + * ( + npa.conj(inputs[0][0]) * inputs[0][3] + - npa.conj(inputs[0][1]) * inputs[0][2] + ) + ) + ) + return flux + + # returns 0,1, or 2 corresponding to x, y, or z normal vectors + # TODO: Handle user-specified normal vectors and cases when 2d + # has a zero-size dimension other than z + def get_normal(self, volume): + # I'll add cylindrical later (since the normal vector gets a little different) + if volume.size.x == 0: return 0 - elif(volume.size.y == 0): + elif volume.size.y == 0: return 1 else: return 2 - diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index 707eb8beb..0a83f0e6e 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -45,6 +45,7 @@ def __init__( self.objective_functions = objective_functions else: self.objective_functions = [objective_functions] + self.objective_arguments = objective_arguments self.f_bank = [] # objective function evaluation history @@ -197,6 +198,7 @@ def forward_run(self): # record objective quantities from user specified monitors self.results_list = [m() for m in self.objective_arguments] + # evaluate objectives self.f0 = [fi(*self.results_list) for fi in self.objective_functions] if len(self.f0) == 1: diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index 5f1935e71..14cb6654f 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -13,7 +13,8 @@ from autograd import tensor_jacobian_product from utils import ApproxComparisonTestCase -MonitorObject = Enum("MonitorObject", "EIGENMODE DFT LDOS POYNTING") +# MonitorObject = Enum("MonitorObject", "EIGENMODE DFT LDOS POYNTING") +MonitorObject = Enum("MonitorObject", "POYNTING") class TestAdjointSolver(ApproxComparisonTestCase): @@ -23,7 +24,7 @@ def setUpClass(cls): cls.silicon = mp.Medium(epsilon=12) cls.sapphire = mp.Medium( - epsilon_diag=(10.225, 10.225, 9.95), + epsilon_cddiag=(10.225, 10.225, 9.95), epsilon_offdiag=(-0.825, -0.55 * np.sqrt(3 / 2), 0.55 * np.sqrt(3 / 2)), ) @@ -169,17 +170,17 @@ def J(mode_mon): def J(mode_mon): return npa.array(mode_mon) - + elif mon_type.name == "POYNTING": obj_list = [ mpa.PoyntingFlux( - sim, - mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0, 1, 0)) + sim, mp.Volume(center=mp.Vector3(1.25), size=mp.Vector3(0, 1, 0)) ) ] - - def J(mode_mon): - return mode_mon + # inputs contains H1, H2, E1, E2, and metadata (in that order) + def J(*inputs): + flux = obj_list[0].compute_flux(*inputs) + return flux opt = mpa.OptimizationProblem( simulation=sim, @@ -342,7 +343,7 @@ def mapping(self, x, filter_radius, eta, beta): projected_field = mpa.tanh_projection(filtered_field, beta, eta) return projected_field.flatten() - + def test_Poynting_Flux(self): print("*** TESTING POYNTING OBJECTIVE ***") @@ -350,12 +351,15 @@ def test_Poynting_Flux(self): for frequencies in [[self.fcen], [1 / 1.58, self.fcen, 1 / 1.53]]: # compute objective value and its gradient for unperturbed design unperturbed_val, unperturbed_grad = self.adjoint_solver( - self.p, MonitorObject.POYNTING, frequencies + self.p, MonitorObject.POYNTING, frequencies, need_gradient=False ) # compute objective value for perturbed design perturbed_val = self.adjoint_solver( - self.p + self.dp, MonitorObject.POYNTING, frequencies, need_gradient=False + self.p + self.dp, + MonitorObject.POYNTING, + frequencies, + need_gradient=False, ) # compare directional derivative @@ -373,7 +377,7 @@ def test_Poynting_Flux(self): tol = 0.07 if mp.is_single_precision() else 0.006 self.assertClose(adj_dd, fnd_dd, epsilon=tol) - + def test_eigenmode(self): print("*** TESTING EIGENMODE OBJECTIVE ***") @@ -434,7 +438,7 @@ def test_DFT_fields(self): tol = 0.07 if mp.is_single_precision() else 0.006 self.assertClose(adj_dd, fnd_dd, epsilon=tol) - + def test_eigenmode(self): print("*** TESTING EIGENMODE OBJECTIVE ***") From 413d154cf20fde61457b92e23db9d901f39b7566 Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Thu, 1 Sep 2022 13:19:30 -0400 Subject: [PATCH 09/10] Still only works in 2D. Fixed cases when some field components are zero (due to symmetry, etc.). Tested with eigenmode and point sources radiating toward line monitors. 3D returns the wrong gradient entirely - I should be able to fix it since it should just be the format of the input. There are a couple cases in 2D simulations where the gradient is somewhat wrong. --- python/adjoint/objective.py | 33 +++++++++++++++++++++++--- python/adjoint/optimization_problem.py | 2 ++ 2 files changed, 32 insertions(+), 3 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 54f94adeb..8f7239cba 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -546,26 +546,52 @@ def register_monitors(self, frequencies): def place_adjoint_source(self, dJ): source = [] + print("This is dJ[0]") + print(dJ[0]) + print(dJ[0].shape) + print("This is metadata's shape:") + print(self.field_component_evaluations[4].shape) for pos, field in enumerate(self.F_fields_list): - source.append(field.place_adjoint_source(np.flipud(dJ[pos].flatten()))[0]) + # Make sure there's a nonzero value in the gradient + # (zero sources don't converge) + # Check is also in prepare_adjoint_run, + # but necessary here too + if np.any(dJ[pos]): + source.append( + field.place_adjoint_source(np.flipud(np.array(dJ[pos]).squeeze()))[ + 0 + ] + ) return source def __call__(self): self.field_component_evaluations = [] + # Get integration weights Meep uses + self.metadata = self.sim.get_array_metadata(vol=self.volume) for field in self.F_fields_list: # Get the dft evaluation from a call to the underlying # FourierFields object field_here = field() + # make sure it isn't a list of scalar zeros equal to the number of + # frequencies (usually caused by symmetries making fields 0) + # fixes the np.array error in the return + # when we give it a "ragged" array + if (np.squeeze(field_here).size) == self._frequencies.size: + print("does the empty array check work") + print(field_here) + field_here = np.array([np.zeros(self.metadata[3].shape)]) self.field_component_evaluations.append(field_here) - # Get integration weights Meep uses - self.metadata = self.sim.get_array_metadata(vol=self.volume) + self.field_component_evaluations.append( np.array([self.metadata[3]]).astype(complex) ) [H1, H2, E1, E2, meta] = self.field_component_evaluations + self._eval = self.field_component_evaluations print("This is meta*E2") print(meta * E2) + print("This is the np array") + print(np.array([H1, H2, E1, E2, meta])) return np.array([H1, H2, E1, E2, meta]) # takes in a 1x5xNxM vector where the size five array corresponds to @@ -573,6 +599,7 @@ def __call__(self): # multiple frequencies will be tested later @staticmethod def compute_flux(*inputs): + # Check if all the inputs are nonzero flux = npa.sum( npa.real( inputs[0][4] diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index 0a83f0e6e..0b06af20b 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -217,6 +217,8 @@ def prepare_adjoint_run(self): for mi, m in enumerate(self.objective_arguments): dJ = jacobian(self.objective_functions[ar], mi)(*self.results_list) # get gradient of objective w.r.t. monitor + print(np.any(dJ)) + print(dJ) if np.any(dJ): self.adjoint_sources[ar] += m.place_adjoint_source( dJ From 1ac0ae1e9d1402473f7e3e30d6d0dbbcc8e20053 Mon Sep 17 00:00:00 2001 From: Alex Kaylor Date: Thu, 15 Sep 2022 15:29:13 -0400 Subject: [PATCH 10/10] Fixed a bug where we take the 0th element the return of FourierFields's place_adjoint_source function. This works in 2D, but it returns 9 separate sources in 3D. The extra sources are handled and collated correctly now. --- python/adjoint/objective.py | 39 +++++++++++++++++++++++++++++-------- 1 file changed, 31 insertions(+), 8 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 8f7239cba..fd5d5fcd5 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -546,23 +546,46 @@ def register_monitors(self, frequencies): def place_adjoint_source(self, dJ): source = [] - print("This is dJ[0]") - print(dJ[0]) - print(dJ[0].shape) + print("This is dJ[0]'s shape:") + print(np.array(dJ[0]).shape) + squeezed_dJ_0 = np.array(dJ[0]).squeeze() + print("This is squeezed dj0 shape:") + print(squeezed_dJ_0.shape) + print("This is metadata's shape:") print(self.field_component_evaluations[4].shape) for pos, field in enumerate(self.F_fields_list): # Make sure there's a nonzero value in the gradient # (zero sources don't converge) # Check is also in prepare_adjoint_run, - # but necessary here too + # but necessary here too since the source is a vector if np.any(dJ[pos]): + reshaped_dJ = np.reshape( + np.array(dJ[pos]).squeeze(), + self.field_component_evaluations[4].shape, + ) + # new_source = field.place_adjoint_source(reshaped_dJ) + new_source = field.place_adjoint_source(np.array(dJ[pos]).squeeze()) + print("This is the new source:") + print(new_source) + print("new soruce's shape") + print(np.array(new_source).shape) source.append( - field.place_adjoint_source(np.flipud(np.array(dJ[pos]).squeeze()))[ - 0 - ] + # field.place_adjoint_source(np.flipud(np.array(dJ[pos]).squeeze()))[ + # 0 + # ] + new_source ) - return source + final_array = np.array(source).flatten() + print("This is the final array shape:") + # print(final_array) + print(np.array(final_array).shape) + # print("This is the final_array with an extra array") + + # test_arr = [final_array] + # print(test_arr) + # print(test_arr.shape) + return final_array.tolist() def __call__(self): self.field_component_evaluations = []