Adjoint near2far with vol_src

python/adjoint/__init__.py

@@ -6,7 +6,7 @@
 
 import meep as mp
 
-from .objective import EigenmodeCoefficient
+from .objective import *
 
 from .basis import BilinearInterpolationBasis
 

python/adjoint/filter_source.py

@@ -11,22 +11,24 @@ def __init__(self,center_frequency,frequencies,frequency_response,dt,time_src=No
         self.center_frequencies = frequencies
 
         # For now, the basis functions cannot overlap much in the frequency domain. Otherwise, the
-        # resulting nodes are wildly large and induce numerical precision errors. We can always 
+        # resulting nodes are wildly large and induce numerical precision errors. We can always
         # produce a safe simulation by forcing the length of each basis function to meet the minimum
         # frequency requirements. This method still minimizes storage requirements.
         self.T = np.max(np.abs(1/np.diff(frequencies)))
         self.N = np.rint(self.T/self.dt)
         self.t = np.arange(0,dt*(self.N),dt)
         self.n = np.arange(self.N)
         f = self.func()
-
+        self.bf = [lambda t, i=i: 0 if t>self.T else (self.nuttall(t,self.center_frequencies)/(self.dt/np.sqrt(2*np.pi)))[i] for i in range(len(self.center_frequencies))]
+        self.time_src_bf = [CustomSource(src_func=bfi,center_frequency=center_frequency,is_integrated=False,end_time=self.T) for bfi in self.bf]
+
         if time_src:
             # get the cutoff of the input signal
             signal_t = np.array([time_src.swigobj.current(ti,dt) for ti in self.t]) # time domain signal
-            signal_dtft = self.dtft(signal_t,self.frequencies)       
+            signal_dtft = self.dtft(signal_t,self.frequencies)
         else:
             signal_dtft = 1
-        
+
         # multiply sampled dft of input signal with filter transfer function
         H = signal_dtft * frequency_response
 
@@ -67,24 +69,25 @@ def nuttall_dtft(self,f,f0):
         return self.cos_window_fd(a,f,f0)
     def dtft(self,y,f):
         return np.matmul(np.exp(1j*2*np.pi*f[:,np.newaxis]*np.arange(y.size)*self.dt), y)*self.dt/np.sqrt(2*np.pi)
-    
+
     def __call__(self,t):
         if t > self.T:
             return 0
         vec = self.nuttall(t,self.center_frequencies) / (self.dt/np.sqrt(2*np.pi)) # compensate for meep dtft
         return np.inner(vec,self.nodes)
-    
+
     def func(self):
-        def _f(t): 
+        def _f(t):
             return self(t)
         return _f
-
+
+
     def estimate_impulse_response(self,H):
         # Use vandermonde matrix to calculate weights of each gaussian. Each window is centered at each frequency point.
         # TODO: come up with a more sophisticated way to choose temporal window size and basis locations
         # that will minimize l2 estimation error and the node weights (since matrix is ill-conditioned)
         vandermonde = self.nuttall_dtft(self.frequencies[:,np.newaxis],self.center_frequencies[np.newaxis,:])
-        nodes = np.matmul(linalg.pinv(vandermonde), H)
+        nodes = np.matmul(linalg.pinv(vandermonde), H.T)
         H_hat = np.matmul(vandermonde, nodes)
-        l2_err = np.sum(np.abs(H-H_hat)**2/np.abs(H)**2)        
+        l2_err = np.sum(np.abs(H-H_hat.T)**2/np.abs(H)**2)
         return nodes, l2_err

python/adjoint/objective.py

@@ -4,6 +4,7 @@
 import numpy as np
 import meep as mp
 from .filter_source import FilteredSource
+from .optimization_problem import atleast_3d, Grid
 
 class ObjectiveQuantitiy(ABC):
     @abstractmethod
@@ -35,15 +36,15 @@ def __init__(self,sim,volume,mode,forward=True,kpoint_func=None,**kwargs):
         self.eval = None
         self.EigenMode_kwargs = kwargs
         return
-    
+
     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.normal_direction = self.monitor.normal_direction
         return self.monitor
-    
+
     def place_adjoint_source(self,dJ):
-        '''Places an equivalent eigenmode monitor facing the opposite direction. Calculates the 
+        '''Places an equivalent eigenmode monitor facing the opposite direction. Calculates the
         correct scaling/time profile.
 
         dJ ........ the user needs to pass the dJ/dMonitor evaluation
@@ -64,31 +65,31 @@ def place_adjoint_source(self,dJ):
                 k0 == direction_scalar * mp.Vector3(z=1)
         else:
             k0 = direction_scalar * self.kpoint_func(self.time_src.frequency,1)
-        
+
         # -------------------------------------- #
-        # Get scaling factor 
+        # Get scaling factor
         # -------------------------------------- #
         # leverage linearity and combine source for multiple frequencies
         if dJ.ndim == 2:
             dJ = np.sum(dJ,axis=1)
-        
+
         # Determine the correct resolution scale factor
         if self.sim.cell_size.y == 0:
             dV = 1/self.sim.resolution
         elif self.sim.cell_size.z == 0:
             dV = 1/self.sim.resolution * 1/self.sim.resolution
         else:
             dV = 1/self.sim.resolution * 1/self.sim.resolution * 1/self.sim.resolution
-        
+
         da_dE = 0.5 * self.cscale # scalar popping out of derivative
         iomega = (1.0 - np.exp(-1j * (2 * np.pi * self.frequencies) * dt)) * (1.0 / dt) # scaled frequency factor with discrete time derivative fix
 
         src = self.time_src
-        
+
         # an ugly way to calcuate the scaled dtft of the forward source
         y = np.array([src.swigobj.current(t,dt) for t in np.arange(0,T,dt)]) # time domain signal
         fwd_dtft = np.matmul(np.exp(1j*2*np.pi*self.frequencies[:,np.newaxis]*np.arange(y.size)*dt), y)*dt/np.sqrt(2*np.pi) # dtft
-        
+
         # we need to compensate for the phase added by the time envelope at our freq of interest
         src_center_dtft = np.matmul(np.exp(1j*2*np.pi*np.array([src.frequency])[:,np.newaxis]*np.arange(y.size)*dt), y)*dt/np.sqrt(2*np.pi)
         adj_src_phase = np.exp(1j*np.angle(src_center_dtft))
@@ -101,7 +102,7 @@ def place_adjoint_source(self,dJ):
             scale = da_dE * dV * dJ * iomega / adj_src_phase
             src = FilteredSource(self.time_src.frequency,self.frequencies,scale,dt) # generate source from broadband response
             amp = 1
-        
+
         # generate source object
         self.source = [mp.EigenModeSource(src,
                     eig_band=self.mode,
@@ -112,7 +113,7 @@ def place_adjoint_source(self,dJ):
                     size=self.volume.size,
                     center=self.volume.center,
                     **self.EigenMode_kwargs)]
-        
+
         return self.source
 
     def __call__(self):
@@ -122,7 +123,7 @@ def __call__(self):
         # Eigenmode data
         direction = mp.NO_DIRECTION if self.kpoint_func else mp.AUTOMATIC
         ob = self.sim.get_eigenmode_coefficients(self.monitor,[self.mode],direction=direction,kpoint_func=self.kpoint_func,**self.EigenMode_kwargs)
-        self.eval = np.squeeze(ob.alpha[:,:,self.forward]) # record eigenmode coefficients for scaling 
+        self.eval = np.squeeze(ob.alpha[:,:,self.forward]) # record eigenmode coefficients for scaling
         self.cscale = ob.cscale # pull scaling factor
 
         return self.eval
@@ -132,4 +133,149 @@ def get_evaluation(self):
         try:
             return self.eval
         except AttributeError:
-            raise RuntimeError("You must first run a forward simulation before resquesting an eigenmode coefficient.")
+            raise RuntimeError("You must first run a forward simulation before resquesting an eigenmode coefficient.")
+
+
+
+class Fourier_Coefficients(ObjectiveQuantitiy):
+    def __init__(self,sim,volume, component):
+
+        self.sim = sim
+        self.volume=volume
+        self.eval = None
+        self.component = component
+        return
+
+    def register_monitors(self,frequencies):
+        self.frequencies = np.asarray(frequencies)
+        self.num_freq = len(self.frequencies)
+        self.monitor = self.sim.add_dft_fields([self.component], self.frequencies, where=self.volume, yee_grid=False)
+
+        return self.monitor
+
+    def place_adjoint_source(self,dJ):
+
+        dt = self.sim.fields.dt
+
+        if self.sim.cell_size.y == 0:
+            dV = 1/self.sim.resolution
+        elif self.sim.cell_size.z == 0:
+            dV = 1/self.sim.resolution * 1/self.sim.resolution
+        else:
+            dV = 1/self.sim.resolution * 1/self.sim.resolution * 1/self.sim.resolution
+
+        self.sources = []
+
+        if self.num_freq == 1:
+            scale = dV * 1j * 2 * np.pi * self.frequencies[0] / self.time_src.fourier_transform(self.frequencies[0])
+            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(self.time_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)):
+                        scale = -dJ_4d[:,xi,yi,zi] * dV * 1j * 2 * np.pi * self.frequencies / np.array([self.time_src.fourier_transform(f) for f in self.frequencies])
+                        src = FilteredSource(self.time_src.frequency,self.frequencies,scale,dt,self.time_src)
+                        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
+
+    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])
+        return self.eval
+
+    def get_evaluation(self):
+
+        try:
+            return self.eval
+        except AttributeError:
+            raise RuntimeError("You must first run a forward simulation.")
+
+
+
+class Far_Coefficients(ObjectiveQuantitiy):
+    def __init__(self,sim,Near2FarRegions, far_pt):
+        self.sim = sim
+        self.Near2FarRegions=Near2FarRegions
+        self.eval = None
+        self.far_pt = far_pt
+        return
+
+    def register_monitors(self,frequencies):
+        self.frequencies = np.asarray(frequencies)
+        self.num_freq = len(self.frequencies)
+        self.monitor = self.sim.add_near2far(self.frequencies, *self.Near2FarRegions, yee_grid=True)
+        return self.monitor
+
+    def place_adjoint_source(self,dJ):
+
+        dt = self.sim.fields.dt
+
+        if self.sim.cell_size.y == 0:
+            dV = 1/self.sim.resolution
+        elif self.sim.cell_size.z == 0:
+            dV = 1/self.sim.resolution * 1/self.sim.resolution
+        else:
+            dV = 1/self.sim.resolution * 1/self.sim.resolution * 1/self.sim.resolution
+
+        self.sources = ['src_vol']
+
+
+        freq_scale = 1j * 2 * np.pi * self.frequencies / np.array([self.time_src.fourier_transform(f) for f in self.frequencies])
+
+        dJ = list(dJ.flatten())
+        dJ_c = mp.ComplexVector(len(dJ))
+        for i in range(len(dJ)):
+            dJ_c[i] = dJ[i]
+
+        #TODO far_pts in 3d or cylindrical, perhaps py_v3_to_vec from simulation.py
+        self.all_nearsrcdata = self.monitor.swigobj.near_sourcedata(mp.vec(self.far_pt.x, self.far_pt.y), dJ_c)
+
+
+        for near_data in self.all_nearsrcdata:
+            cur_comp = near_data.near_fd_comp
+            amp_arr = near_data.amp_arr.reshape(self.num_freq, -1)
+
+
+            if cur_comp in [mp.Ex, mp.Ey, mp.Ez]:
+                scale = -freq_scale * amp_arr
+            else:
+                scale = freq_scale * amp_arr
+
+            if self.num_freq == 1:
+                ##TODO
+                self.sources += [mp.Source(self.time_src, component=cur_comp, amplitude=scale[0], center=near_pt)]
+            else:
+                src = FilteredSource(self.time_src.frequency,self.frequencies,scale,dt,self.time_src)
+                self.sources+=[(src.nodes, src.time_src_bf, near_data)]
+
+
+
+
+        return self.sources
+
+    def __call__(self):
+        self.time_src = self.sim.sources[0].src
+        self.eval = np.array(self.sim.get_farfield(self.monitor, self.far_pt))
+        self.eval = self.eval.reshape((self.num_freq, 6))
+        return self.eval
+
+    def get_evaluation(self):
+        try:
+            return self.eval
+        except AttributeError:
+            raise RuntimeError("You must first run a forward simulation.")