WIP: Adjoint poynting flux

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
 
 
 

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
 
@@ -70,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
@@ -124,38 +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)
 
+        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:
-            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]))]
+            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:
-            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
+            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):
@@ -213,6 +222,77 @@ 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):        
+
+        # 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
+
+    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 +301,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:

python/meep.i

@@ -1479,7 +1479,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
 %template(ComplexVector) std::vector<std::complex<double> >;
 
 std::vector<struct meep::sourcedata> meep::dft_near2far::near_sourcedata(const meep::vec &x, std::complex<double>* dJ);
-
+std::vector<struct meep::sourcedata> meep::dft_flux::flux_sourcedata(std::complex<double>* dJ);
 void meep::fields::add_srcdata(struct meep::sourcedata cur_data, meep::src_time *src, size_t n, std::complex<double>* amp_arr);
 
 struct vector3 {

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)
 

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):
@@ -587,8 +593,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
+

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;