From f228b68bb31d43fef95209c26f047d1d0ed0cbdc Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Fri, 9 Oct 2020 10:47:56 -0400 Subject: [PATCH 1/4] fix some adjoint filter bugs --- python/adjoint/filters.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 2d632162d..ef79be379 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -5,6 +5,7 @@ import numpy as np from autograd import numpy as npa import meep as mp +from scipy import special def _centered(arr, newshape): '''Helper function that reformats the padded array of the fft filter operation. @@ -581,7 +582,7 @@ def tanh_projection(x,beta,eta): return (npa.tanh(beta*eta) + npa.tanh(beta*(x-eta))) / (npa.tanh(beta*eta) + npa.tanh(beta*(1-eta))) -def heaviside_projection(x, beta): +def heaviside_projection(x, beta, eta): '''Projection filter that thresholds the input parameters between 0 and 1. Parameters @@ -880,7 +881,7 @@ def gray_indicator(x): [1] Lazarov, B. S., Wang, F., & Sigmund, O. (2016). Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. ''' - return np.mean(4 * x.flatten()) * (1-x.flatten()) * 100 + return npa.mean(4 * x.flatten() * (1-x.flatten())) * 100 From c44471abd517266ca4ab5e900c1bc9321647d214 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Fri, 9 Oct 2020 17:21:07 -0400 Subject: [PATCH 2/4] framework --- python/adjoint/objective.py | 66 ++++++++++++++++++++++++++++++++++++- python/meep.i | 2 +- python/source.py | 5 ++- src/dft.cpp | 65 ++++++++++++++++++++++++++++++++++++ src/meep.hpp | 18 +++++----- 5 files changed, 144 insertions(+), 12 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 7d62511b1..d3ac2b0ff 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -6,6 +6,9 @@ from .filter_source import FilteredSource from .optimization_problem import atleast_3d, Grid +# Transverse component definitions needed for flux adjoint calculations +EH_components = [[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 ObjectiveQuantitiy(ABC): @abstractmethod def __init__(self): @@ -39,7 +42,7 @@ def __init__(self,sim,volume,mode,forward=True,kpoint_func=None,**kwargs): def register_monitors(self,frequencies): self.frequencies = np.asarray(frequencies) - self.monitor = self.sim.add_mode_monitor(frequencies,mp.FluxRegion(center=self.volume.center,size=self.volume.size),yee_grid=True) + self.monitor = self.sim.add_mode_monitor(frequencies,mp.ModeRegion(center=self.volume.center,size=self.volume.size),yee_grid=True) self.normal_direction = self.monitor.normal_direction return self.monitor @@ -213,6 +216,61 @@ def get_evaluation(self): except AttributeError: raise RuntimeError("You must first run a forward simulation.") +class PoyntingFlux(ObjectiveQuantitiy): + def __init__(self,sim,volume): + self.sim = sim + self.volume=volume + self.eval = None + self.normal_direction = None + return + + def register_monitors(self,frequencies): + self.frequencies = np.asarray(frequencies) + self.num_freq = len(self.frequencies) + self.monitor = self.sim.add_flux(self.frequencies, mp.FluxRegion(center=self.volume.center,size=self.volume.size)) + self.normal_direction = self.monitor.normal_direction + return self.monitor + + def place_adjoint_source(self,dJ): + ''' + We'll add one volume source for every chunk at every frequency. + All of the scaling will occur inside of the `flux_sourcedata()` + function. So long as we cache the results of our temporal + basis function, this should be rather computationally cheap. + ''' + dt = self.sim.fields.dt # the timestep size from sim.fields.dt of the forward sim + self.sources = [] + dJ = np.atleast_1d(dJ) + dJ = np.sum(dJ,axis=1) + + # Pull a list of amp array data for each chunk on the current proc + self.all_fluxsrcdata = self.monitor.swigobj.flux_sourcedata(dJ) + # Loop over chunk list and create a volume source for each chunk + for flux_data in self.all_fluxsrcdata: + amp_arr = np.array(flux_data.amp_arr).reshape(-1, self.num_freq) + scale = amp_arr * adj_src_scale(self, dt) + if self.num_freq == 1: + self.sources += [mp.IndexedSource(self.time_src, flux_data, scale[:,0],center=self.volume.center,size=self.volume.size)] + else: + src = FilteredSource(self.time_src.frequency,self.frequencies,scale,dt) + (num_basis, num_pts) = src.nodes.shape + # For each frequency (i.e. basis function) add a volume source. + for basis_i in range(num_basis): + self.sources += [mp.IndexedSource(src.time_src_bf[basis_i], flux_data, src.nodes[basis_i],center=self.volume.center,size=self.volume.size)] + + return self.sources + + def __call__(self): + self.time_src = self.sim.sources[0].src + self.eval = np.array(mp.get_fluxes(self.monitor)) + + return self.eval + + def get_evaluation(self): + try: + return self.eval + except AttributeError: + raise RuntimeError("You must first run a forward simulation.") def adj_src_scale(obj_quantity, dt, include_resolution=True): # -------------------------------------- # @@ -221,6 +279,12 @@ def adj_src_scale(obj_quantity, dt, include_resolution=True): # leverage linearity and combine source for multiple frequencies T = obj_quantity.sim.meep_time() + ''' + Integral-like adjoints (e.g. poynting flux, mode overlaps, etc) + have a dV scale factor that needs to be included in the adjoint + source. Other quantities (e.g. Near2Far, DFT fields, etc.) don't + need the scale factor since no integral is involved. + ''' if not include_resolution: dV = 1 elif obj_quantity.sim.cell_size.y == 0: diff --git a/python/meep.i b/python/meep.i index 61354bb95..3282cc703 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1479,7 +1479,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj %template(ComplexVector) std::vector >; std::vector meep::dft_near2far::near_sourcedata(const meep::vec &x, std::complex* dJ); - +std::vector meep::dft_flux::flux_sourcedata(std::complex* dJ); void meep::fields::add_srcdata(struct meep::sourcedata cur_data, meep::src_time *src, size_t n, std::complex* amp_arr); struct vector3 { diff --git a/python/source.py b/python/source.py index 121ac00a1..25f10c4c0 100644 --- a/python/source.py +++ b/python/source.py @@ -587,8 +587,11 @@ class IndexedSource(Source): """ created a source object using (SWIG-wrapped mp::srcdata*) srcdata. """ - def __init__(self, src, srcdata, amp_arr): + def __init__(self, src, srcdata, amp_arr,vol=None, center=None, size=None): self.src = src self.num_pts = len(amp_arr) self.srcdata = srcdata self.amp_arr = amp_arr + self.center=center + self.size=size + diff --git a/src/dft.cpp b/src/dft.cpp index c3e392fa8..9a846dbb2 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include "meep.hpp" #include "meep_internals.hpp" @@ -1249,4 +1250,68 @@ void fields::get_mode_mode_overlap(void *mode1_data, void *mode2_data, dft_flux get_overlap(mode1_data, mode2_data, flux, 0, overlaps); } +//Modified from dft_near2far::near_sourcedata +std::vector dft_flux::flux_sourcedata(std::complex *dJ) { + + const size_t Nfreq = freq.size(); + std::vector temp; + + // Determine which components we care about + component cE_t[2] = {Ex, Ey}, cH_t[2] = {Hy, Hx}; + switch (normal_direction) { + case X: cE_t[0] = Ey, cE_t[1] = Ez, cH_t[0] = Hz, cH_t[1] = Hy; break; + case Y: cE_t[0] = Ez, cE_t[1] = Ex, cH_t[0] = Hx, cH_t[1] = Hz; break; + case R: cE_t[0] = Ep, cE_t[1] = Ez, cH_t[0] = Hz, cH_t[1] = Hp; break; + case P: cE_t[0] = Ez, cE_t[1] = Er, cH_t[0] = Hr, cH_t[1] = Hz; break; + case Z: + if (cE==Er) // check the stored cE, which only stores the first component + cE_t[0] = Er, cE_t[1] = Ep, cH_t[0] = Hp, cH_t[1] = Hr; + else + cE_t[0] = Ex, cE_t[1] = Ey, cH_t[0] = Hy, cH_t[1] = Hx; + break; + default: abort("invalid flux component!"); + } + component c0[4] = {cE_t[0], cE_t[1], cH_t[0], cH_t[1]}; + component c_adj[4] = {cH_t[1], cH_t[0], cE_t[1], cE_t[0]}; + + // Loop over chunkloops and chunks + dft_chunk *chunklists[2]; + chunklists[0] = E; + chunklists[1] = H; + for (int ncl = 0; ncl < 2; ncl++) { + int chunk_idx = 0; + for (dft_chunk *chunk = chunklists[ncl]; chunk; chunk = chunk->next_in_dft) { + int ci = ncl * 2 + chunk_idx; + component currentComp = c0[ci]; + component adjointComp = c_adj[ci]; + + assert(Nfreq == chunk->omega.size()); + std::vector idx_arr; + std::vector > amp_arr; + + vec rshift(chunk->shift * (0.5 * chunk->fc->gv.inva)); + size_t idx_dft = 0; + sourcedata temp_struct = {adjointComp, idx_arr, chunk->fc->chunk_idx, amp_arr}; + LOOP_OVER_IVECS(chunk->fc->gv, chunk->is, chunk->ie, idx) { + IVEC_LOOP_LOC(chunk->fc->gv, x0); + x0 = chunk->S.transform(x0, chunk->sn) + rshift; + vec xs(x0); + double w = IVEC_LOOP_WEIGHT(chunk->s0, chunk->s1, chunk->e0, chunk->e1, chunk->dV0 + chunk->dV1 * loop_i2); + + temp_struct.idx_arr.push_back(idx); + for (size_t i = 0; i < Nfreq; ++i) { + cdouble dft_val = chunk->dft[chunk->omega.size() * (idx_dft++) + Nfreq]/w; // chunk->stored_weight; + dft_val = std::conj(dft_val); + // principle of equivalence + if (is_electric(currentComp)) dft_val *= -1; + + temp_struct.amp_arr.push_back(dft_val); + } + } + temp.push_back(temp_struct); + } + } + return temp; +} + } // namespace meep diff --git a/src/meep.hpp b/src/meep.hpp index 59cfc31fb..ba0b00e1f 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -995,6 +995,13 @@ class monitor_point { double fmin, double fmax, int maxbands); }; +struct sourcedata{ + component near_fd_comp; + std::vector idx_arr; + int fc_idx; + std::vector > amp_arr; +}; + // dft.cpp // this should normally only be created with fields::add_dft class dft_chunk { @@ -1125,6 +1132,8 @@ class dft_flux { volume where; direction normal_direction; bool use_symmetry; + + std::vector flux_sourcedata(std::complex* dJ); }; // dft.cpp (normally created with fields::add_dft_energy) @@ -1194,15 +1203,6 @@ class dft_force { volume where; }; - -struct sourcedata{ - component near_fd_comp; - std::vector idx_arr; - int fc_idx; - std::vector > amp_arr; -}; - - // near2far.cpp (normally created with fields::add_dft_near2far) class dft_near2far { public: From 59ca50e11fc77894fd118a8f8a0347deaa18fb61 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Mon, 12 Oct 2020 20:54:19 -0400 Subject: [PATCH 3/4] fixed EM crash and started fourier coefficients --- python/adjoint/objective.py | 40 +++++++++++++---------------- src/dft.cpp | 50 ++++++++++++++++++++++++++++++++++++- src/meep.hpp | 3 +++ 3 files changed, 69 insertions(+), 24 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index d3ac2b0ff..579addb7c 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -73,6 +73,7 @@ def place_adjoint_source(self,dJ): if self.frequencies.size == 1: # Single frequency simulations. We need to drive it with a time profile. amp = da_dE * dJ * scale # final scale factor + src = self.time_src else: # multi frequency simulations scale = da_dE * dJ * scale @@ -130,30 +131,23 @@ def place_adjoint_source(self,dJ): dt = self.sim.fields.dt # the timestep size from sim.fields.dt of the forward sim self.sources = [] scale = adj_src_scale(self, dt) + # Pull a list of amp array data for each chunk on the current proc + self.all_dftsrcdata = self.monitor.swigobj.dft_sourcedata(dJ) + # Loop over chunk list and create a volume source for each chunk + for dft_data in self.all_dftsrcdata: + amp_arr = np.array(dft_data.amp_arr).reshape(-1, self.num_freq) + scale = amp_arr.flatten()*adj_src_scale(self, dt).flatten()*dJ.flatten() + print(scale) + if self.num_freq == 1: + self.sources += [mp.IndexedSource(self.time_src, dft_data, scale,center=self.volume.center,size=self.volume.size)] + else: + src = FilteredSource(self.time_src.frequency,self.frequencies,scale,dt) + (num_basis, num_pts) = src.nodes.shape + # For each frequency (i.e. basis function) add a volume source. + for basis_i in range(num_basis): + self.sources += [mp.IndexedSource(src.time_src_bf[basis_i], flux_data, src.nodes[basis_i],center=self.volume.center,size=self.volume.size)] - if self.num_freq == 1: - amp = -atleast_3d(dJ[0]) * scale - if self.component in [mp.Hx, mp.Hy, mp.Hz]: - amp = -amp - for zi in range(len(self.dg.z)): - for yi in range(len(self.dg.y)): - for xi in range(len(self.dg.x)): - if amp[xi, yi, zi] != 0: - self.sources += [mp.Source(src, component=self.component, amplitude=amp[xi, yi, zi], - center=mp.Vector3(self.dg.x[xi], self.dg.y[yi], self.dg.z[zi]))] - else: - dJ_4d = np.array([atleast_3d(dJ[f]) for f in range(self.num_freq)]) - if self.component in [mp.Hx, mp.Hy, mp.Hz]: - dJ_4d = -dJ_4d - for zi in range(len(self.dg.z)): - for yi in range(len(self.dg.y)): - for xi in range(len(self.dg.x)): - final_scale = -dJ_4d[:,xi,yi,zi] * scale - src = FilteredSource(self.time_src.frequency,self.frequencies,final_scale,dt) - self.sources += [mp.Source(src, component=self.component, amplitude=1, - center=mp.Vector3(self.dg.x[xi], self.dg.y[yi], self.dg.z[zi]))] - - return self.sources + return self.sources def __call__(self): self.time_src = self.sim.sources[0].src diff --git a/src/dft.cpp b/src/dft.cpp index 9a846dbb2..ecd152435 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -1087,7 +1087,8 @@ cdouble *fields::get_dft_array(dft_fields fdft, component c, int num_freq, int * cdouble *array; direction dirs[3]; process_dft_component(chunklists, 1, num_freq, c, 0, &array, rank, dims, dirs); - return collapse_array(array, rank, dims, dirs, fdft.where); + return collapse_array(array, rank, dims, dirs, fdft.where);; + } /***************************************************************/ @@ -1314,4 +1315,51 @@ std::vector dft_flux::flux_sourcedata(std::complex *d return temp; } +//Modified from dft_near2far::near_sourcedata +std::vector dft_fields::dft_sourcedata(std::complex *dJ) { + + const size_t Nfreq = freq.size(); + std::vector temp; + + // Loop over chunkloops and chunks + int c_i = 0; + for (dft_chunk *f = chunks; f; f = f->next_in_dft) { + component currentComp = f->c; + component adjointComp = currentComp; // for simple dft fields, adjoint comp is same as forward comp + + assert(Nfreq == f->omega.size()); + std::vector idx_arr; + std::vector > amp_arr; + + vec rshift(f->shift * (0.5 * f->fc->gv.inva)); + size_t idx_dft = 0; + sourcedata temp_struct = {adjointComp, idx_arr, f->fc->chunk_idx, amp_arr}; + LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) { + + //Calculate the interpolation weights + IVEC_LOOP_LOC(f->fc->gv, x0); + x0 = f->S.transform(x0, f->sn) + rshift; + vec xs(x0); + double w = IVEC_LOOP_WEIGHT(f->s0, f->s1, f->e0, f->e1, f->dV0 + f->dV1 * loop_i2); + + temp_struct.idx_arr.push_back(idx); // store index of current point + + // iterate through frequency points + for (size_t i = 0; i < Nfreq; ++i) { + + //cdouble dft_val = f->dft[f->omega.size() * (idx_dft++) + Nfreq]; + //dft_val = std::conj(dft_val); + cdouble dft_val = 1; + // principle of equivalence + if (is_electric(currentComp)) dft_val *= -1; + + temp_struct.amp_arr.push_back(dft_val); + } + } + temp.push_back(temp_struct); + c_i++; + } + return temp; +} + } // namespace meep diff --git a/src/meep.hpp b/src/meep.hpp index ba0b00e1f..a67d11469 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1304,6 +1304,9 @@ class dft_fields { std::vector freq; dft_chunk *chunks; volume where; + + // compute dft adjoint sources + std::vector dft_sourcedata(std::complex* dJ); }; enum in_or_out { Incoming = 0, Outgoing }; From 0f134f3b41dbc56113cd7338565276ce613cc94c Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 20 Oct 2020 16:27:04 -0400 Subject: [PATCH 4/4] switch to amp_data source approach --- python/adjoint/objective.py | 120 ++++++++++++++++++++++-------------- python/simulation.py | 16 ++++- python/source.py | 8 ++- src/array_slice.cpp | 4 +- src/meep.hpp | 2 +- src/sources.cpp | 46 +++++++++++--- 6 files changed, 135 insertions(+), 61 deletions(-) diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 579addb7c..ca86d3fde 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -128,31 +128,43 @@ def register_monitors(self,frequencies): return self.monitor def place_adjoint_source(self,dJ): + # Correctly format the dJ matrix + print(dJ.shape) + dJ_shape = np.array(self.weights.shape) + dJ_shape[-1] = self.num_freq + dJ = np.ascontiguousarray((dJ).reshape(dJ_shape)) + dt = self.sim.fields.dt # the timestep size from sim.fields.dt of the forward sim - self.sources = [] - scale = adj_src_scale(self, dt) - # Pull a list of amp array data for each chunk on the current proc - self.all_dftsrcdata = self.monitor.swigobj.dft_sourcedata(dJ) - # Loop over chunk list and create a volume source for each chunk - for dft_data in self.all_dftsrcdata: - amp_arr = np.array(dft_data.amp_arr).reshape(-1, self.num_freq) - scale = amp_arr.flatten()*adj_src_scale(self, dt).flatten()*dJ.flatten() - print(scale) - if self.num_freq == 1: - self.sources += [mp.IndexedSource(self.time_src, dft_data, scale,center=self.volume.center,size=self.volume.size)] - else: - src = FilteredSource(self.time_src.frequency,self.frequencies,scale,dt) - (num_basis, num_pts) = src.nodes.shape - # For each frequency (i.e. basis function) add a volume source. - for basis_i in range(num_basis): - self.sources += [mp.IndexedSource(src.time_src_bf[basis_i], flux_data, src.nodes[basis_i],center=self.volume.center,size=self.volume.size)] - return self.sources + mon_dv = self.sim.resolution ** self.volume.get_nonzero_dims() + time_scale = adj_src_scale(self, dt) + amp = dJ*time_scale*mon_dv + self.sources = [] + if self.component in [mp.Ex, mp.Ey, mp.Ez]: + amp = -amp + if self.num_freq == 1: + self.sources += [mp.Source(self.time_src,component=self.component,amp_data=amp[:,:,:,0], + center=self.volume.center,size=self.volume.size,amp_data_use_grid=True)] + else: + src = FilteredSource(self.time_src.frequency,self.frequencies,amp.reshape(-1, *amp.shape[-1:]),dt) + (num_basis, num_pts) = src.nodes.shape + fit_data = src.nodes + fit_data = (fit_data.T)[:,np.newaxis,np.newaxis,:]#fit_data.reshape(dJ_shape) + for basis_i in range(num_basis): + self.sources += [mp.Source(src.time_src_bf[basis_i],component=self.component,amp_data=np.ascontiguousarray(fit_data[:,:,:,basis_i]), + center=self.volume.center,size=self.volume.size,amp_data_use_grid=True)] + return self.sources def __call__(self): self.time_src = self.sim.sources[0].src - self.dg = Grid(*self.sim.get_array_metadata(dft_cell=self.monitor)) self.eval = np.array([self.sim.get_dft_array(self.monitor, self.component, i) for i in range(self.num_freq)]) #Shape = (num_freq, [pts]) + + # Calculate, shape, and store weights + self.dg = Grid(*self.sim.get_array_metadata(dft_cell=self.monitor)) + self.expected_dims = np.array([len(self.dg.x),len(self.dg.y),len(self.dg.z)]) + self.weights = self.dg.w.reshape(self.expected_dims) + self.weights = self.weights[...,np.newaxis] + return self.eval def get_evaluation(self): @@ -225,38 +237,54 @@ def register_monitors(self,frequencies): self.normal_direction = self.monitor.normal_direction return self.monitor - def place_adjoint_source(self,dJ): - ''' - We'll add one volume source for every chunk at every frequency. - All of the scaling will occur inside of the `flux_sourcedata()` - function. So long as we cache the results of our temporal - basis function, this should be rather computationally cheap. - ''' - dt = self.sim.fields.dt # the timestep size from sim.fields.dt of the forward sim - self.sources = [] - dJ = np.atleast_1d(dJ) - dJ = np.sum(dJ,axis=1) + def place_adjoint_source(self,dJ): - # Pull a list of amp array data for each chunk on the current proc - self.all_fluxsrcdata = self.monitor.swigobj.flux_sourcedata(dJ) - # Loop over chunk list and create a volume source for each chunk - for flux_data in self.all_fluxsrcdata: - amp_arr = np.array(flux_data.amp_arr).reshape(-1, self.num_freq) - scale = amp_arr * adj_src_scale(self, dt) - if self.num_freq == 1: - self.sources += [mp.IndexedSource(self.time_src, flux_data, scale[:,0],center=self.volume.center,size=self.volume.size)] - else: - src = FilteredSource(self.time_src.frequency,self.frequencies,scale,dt) - (num_basis, num_pts) = src.nodes.shape - # For each frequency (i.e. basis function) add a volume source. - for basis_i in range(num_basis): - self.sources += [mp.IndexedSource(src.time_src_bf[basis_i], flux_data, src.nodes[basis_i],center=self.volume.center,size=self.volume.size)] - - return self.sources + # Format jacobian (and use linearity) + if dJ.ndim == 2: + dJ = np.sum(dJ,axis=1) + dJ = dJ.reshape(1,1,1,self.num_freq,1) + + # Adjust for time scaling + dt = self.sim.fields.dt + time_scale = adj_src_scale(self, dt).reshape(1,1,1,self.num_freq,1) + + # final source amplitude as a function of position + amp = dJ * time_scale * np.conj(self.m_EH) * self.weights * (self.sim.resolution**self.volume.get_nonzero_dims()) + + self.sources = [] + for ic,c in enumerate(EH_components[self.normal_direction]): + # principle of equivalence + amp_scale = -1 if c in [mp.Hx,mp.Hy,mp.Hz] else 1 + + if np.any(amp[:,:,:,:,ic]): + # single frequency case + if self.num_freq == 1: + self.sources += [mp.Source(self.time_src,component=EH_components[self.normal_direction][3-ic],amp_data=np.ascontiguousarray(amp[:,:,:,0,ic]), + center=self.volume.center,size=self.volume.size,amp_data_use_grid=True,amplitude=amp_scale)] + # multi frequency case + else: + src = FilteredSource(self.time_src.frequency,self.frequencies,amp[...,ic].reshape(-1, self.num_freq),dt) # fit data to basis functions + fit_data = (src.nodes.T).reshape(amp[...,ic].shape) # format new amplitudes + for basis_i in range(self.num_freq): + self.sources += [mp.Source(src.time_src_bf[basis_i],EH_components[self.normal_direction][3-ic],amp_data=np.ascontiguousarray(fit_data[:,:,:,basis_i]), + center=self.volume.center,size=self.volume.size,amp_data_use_grid=True,amplitude=amp_scale)] + return self.sources def __call__(self): self.time_src = self.sim.sources[0].src + + # Evaluate poynting flux self.eval = np.array(mp.get_fluxes(self.monitor)) + + # Calculate, shape, and store weights + self.dg = Grid(*self.sim.get_array_metadata(dft_cell=self.monitor)) + self.weights = self.dg.w.reshape((len(self.dg.x),len(self.dg.y),len(self.dg.z),1,1)) + + # Store fields at monitor + self.m_EH = np.zeros((len(self.dg.x), len(self.dg.y), len(self.dg.z), self.num_freq, 4), dtype=complex) + for f in range(self.num_freq): + for ic, c in enumerate(EH_components[self.normal_direction]): + self.m_EH[..., f, ic] = self.sim.get_dft_array(self.monitor, c, f).reshape(len(self.dg.x), len(self.dg.y), len(self.dg.z)) return self.eval diff --git a/python/simulation.py b/python/simulation.py index eda662dd7..02d5ff2a6 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -434,6 +434,20 @@ def pt_in_volume(self,pt): return True else: return False + + def get_nonzero_dims(self): + if ((self.size.x != 0) and (self.size.y != 0) and (self.size.z != 0)): + return 3 + elif (((self.size.x == 0) and (self.size.y != 0) and (self.size.z != 0)) or + ((self.size.x != 0) and (self.size.y == 0) and (self.size.z != 0)) or + ((self.size.x != 0) and (self.size.y != 0) and (self.size.z == 0))): + return 2 + elif (((self.size.x == 0) and (self.size.y == 0) and (self.size.z != 0)) or + ((self.size.x != 0) and (self.size.y == 0) and (self.size.z == 0)) or + ((self.size.x == 0) and (self.size.y != 0) and (self.size.z == 0))): + return 1 + else: + return 0 class FluxRegion(object): @@ -2372,7 +2386,7 @@ def add_source(self, src): elif src.amp_func: add_vol_src(src.amp_func, src.amplitude * 1.0) elif src.amp_data is not None: - add_vol_src(src.amp_data, src.amplitude * 1.0) + add_vol_src(src.amp_data, src.amplitude * 1.0, src.amp_data_use_grid) else: add_vol_src(src.amplitude * 1.0) diff --git a/python/source.py b/python/source.py index 25f10c4c0..e73b8e97c 100644 --- a/python/source.py +++ b/python/source.py @@ -37,7 +37,7 @@ class Source(object): Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). """ def __init__(self, src, component, center=None, volume=None, size=Vector3(), amplitude=1.0, amp_func=None, - amp_func_file='', amp_data=None): + amp_func_file='', amp_data=None,amp_data_use_grid=False): """ Construct a `Source`. @@ -83,6 +83,11 @@ def __init__(self, src, component, center=None, volume=None, size=Vector3(), amp For a 2d simulation, just pass 1 for the third dimension, e.g., `arr = np.zeros((N, M, 1), dtype=np.complex128)`. Defaults to `None`. + + **`amp_data_use_grid` [`bool`]** — Default: False. If either `amp_data` or `amp_func_file` + are used, specify whether to use the typical volume edges when interpolating + (the default behavior) or whether to use the grid returned from `get_array_metadata` + (i.e. when using this feature for adjoint calculations). + As described in Section 4.2 ("Incident Fields and Equivalent Currents") in [Chapter 4](http://arxiv.org/abs/arXiv:1301.5366) ("Electromagnetic Wave Source Conditions") of the book [Advances in FDTD Computational Electrodynamics: @@ -113,6 +118,7 @@ def __init__(self, src, component, center=None, volume=None, size=Vector3(), amp self.amp_func = amp_func self.amp_func_file = amp_func_file self.amp_data = amp_data + self.amp_data_use_grid = amp_data_use_grid class SourceTime(object): diff --git a/src/array_slice.cpp b/src/array_slice.cpp index 0c16521b3..4182267c5 100644 --- a/src/array_slice.cpp +++ b/src/array_slice.cpp @@ -452,8 +452,8 @@ int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], dire LOOP_OVER_DIRECTIONS(gv.dim, d) { if (rank >= 3) abort("too many dimensions in array_slice"); size_t n = (data->max_corner.in_direction(d) - data->min_corner.in_direction(d)) / 2 + 1; - if (where.in_direction(d) == 0.0 && collapse_empty_dimensions) n = 1; - if (n > 1) { + if (where.in_direction(d) == 0.0) n = 1; + if (!collapse_empty_dimensions) { data->ds[rank] = d; dims[rank++] = n; slice_size *= n; diff --git a/src/meep.hpp b/src/meep.hpp index a67d11469..c0b90bb93 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1728,7 +1728,7 @@ class fields { void add_volume_source(component c, const src_time &src, const volume &where_, std::complex *arr, size_t dim1, size_t dim2, size_t dim3, - std::complex amp); + std::complex amp, bool use_grid=false); void add_volume_source(component c, const src_time &src, const volume &where_, const char *filename, const char *dataset, std::complex amp); void add_volume_source(component c, const src_time &src, const volume &, diff --git a/src/sources.cpp b/src/sources.cpp index 148121c24..c59aa6859 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -325,7 +325,8 @@ void fields::add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, st static realnum *amp_func_data_re = NULL; static realnum *amp_func_data_im = NULL; static const volume *amp_func_vol = NULL; -static size_t amp_file_dims[3]; +static realnum ampl_func_dims[3]; // effective volume dimensions +static size_t amp_file_dims[3]; // array dimensions complex amp_file_func(const vec &p) { double x_size = 0, y_size = 0, z_size = 0; @@ -333,17 +334,17 @@ complex amp_file_func(const vec &p) { switch (amp_func_vol->dim) { case D1: z_size = amp_func_vol->in_direction(Z); break; case D2: - x_size = amp_func_vol->in_direction(X); - y_size = amp_func_vol->in_direction(Y); + x_size = ampl_func_dims[0]; + y_size = ampl_func_dims[1]; break; case D3: - x_size = amp_func_vol->in_direction(X); - y_size = amp_func_vol->in_direction(Y); - z_size = amp_func_vol->in_direction(Z); + x_size = ampl_func_dims[0]; + y_size = ampl_func_dims[1]; + z_size = ampl_func_dims[2]; break; case Dcyl: - x_size = amp_func_vol->in_direction(X); - z_size = amp_func_vol->in_direction(Z); + x_size = ampl_func_dims[0]; + z_size = ampl_func_dims[2]; break; } @@ -361,7 +362,7 @@ complex amp_file_func(const vec &p) { void fields::add_volume_source(component c, const src_time &src, const volume &where_, complex *arr, size_t dim1, size_t dim2, size_t dim3, - complex amp) { + complex amp, bool use_grid) { amp_func_vol = &where_; @@ -377,7 +378,32 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w amp_func_data_re[i] = real(arr[i]); amp_func_data_im[i] = imag(arr[i]); } - + + /* The default amp_array interpolation behavior + is to use the actual borders of the volume object + as the edges of our interpolation cell. + + In some cases,(e.g. adjoint calculations) it's important + to define our interpolation cell by what's returned by + get_dft_array(). This often corresponds to a volume + that is padded by an extra pixel on each side. We can + consistenly predict the size of this volume using + get_array_slice_dimensions().*/ + if (use_grid){ + size_t dims[3]; + direction dirs[3]; + vec min_max_loc[2]; // extremal points in subgrid + int rank = get_array_slice_dimensions(*amp_func_vol, dims, dirs, false /*collapse_empty_dimensions*/, + true /*snap_empty_dimensions*/, + min_max_loc); + ampl_func_dims[0] = min_max_loc[1].x() - min_max_loc[0].x(); + ampl_func_dims[1] = min_max_loc[1].y() - min_max_loc[0].y(); + ampl_func_dims[2] = min_max_loc[1].z() - min_max_loc[0].z(); + } else{ + ampl_func_dims[0] = amp_func_vol->in_direction(X); + ampl_func_dims[1] = amp_func_vol->in_direction(Y); + ampl_func_dims[2] = amp_func_vol->in_direction(Z); + } add_volume_source(c, src, where_, amp_file_func, amp); delete[] amp_func_data_re;