Add support for load_minus_flux to objective functions of adjoint solver

python/adjoint/objective.py

@@ -4,7 +4,7 @@
 from typing import Callable, List, Optional
 
 import numpy as np
-from meep.simulation import py_v3_to_vec
+from meep.simulation import py_v3_to_vec, FluxData, NearToFarData
 
 import meep as mp
 
@@ -158,6 +158,11 @@ class EigenmodeCoefficient(ObjectiveQuantity):
         kpoint_func_overlap_idx: the index of the mode coefficient to return when
           specifying `kpoint_func`. When specified, this overrides the effect of
           `forward` and should have a value of either 0 or 1.
+        subtracted_dft_fields: the DFT fields obtained using `get_flux_data` from
+          a previous normalization run. This is subtracted from the DFT fields
+          of this mode monitor in order to improve the accuracy of the
+          reflectance measurement (i.e., the $S_{11}$ scattering parameter).
+          Default is None.
     """
 
     def __init__(
@@ -169,6 +174,7 @@ def __init__(
         kpoint_func: Optional[Callable] = None,
         kpoint_func_overlap_idx: Optional[int] = 0,
         decimation_factor: Optional[int] = 0,
+        subtracted_dft_fields: Optional[FluxData] = None,
         **kwargs
     ):
         """
@@ -189,6 +195,7 @@ def __init__(
         self._monitor = None
         self._cscale = None
         self.decimation_factor = decimation_factor
+        self.subtracted_dft_fields = subtracted_dft_fields
 
     def register_monitors(self, frequencies):
         self._frequencies = np.asarray(frequencies)
@@ -198,6 +205,11 @@ def register_monitors(self, frequencies):
             yee_grid=True,
             decimation_factor=self.decimation_factor,
         )
+        if self.subtracted_dft_fields is not None:
+            self.sim.load_minus_flux_data(
+                self._monitor,
+                self.subtracted_dft_fields,
+            )
         return self._monitor
 
     def place_adjoint_source(self, dJ):
@@ -279,13 +291,15 @@ def __init__(
         component: List[int],
         yee_grid: Optional[bool] = False,
         decimation_factor: Optional[int] = 0,
+        subtracted_dft_fields: Optional[FluxData] = None,
     ):
         """ """
         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.subtracted_dft_fields = subtracted_dft_fields
 
     def register_monitors(self, frequencies):
         self._frequencies = np.asarray(frequencies)
@@ -296,6 +310,11 @@ def register_monitors(self, frequencies):
             yee_grid=self.yee_grid,
             decimation_factor=self.decimation_factor,
         )
+        if self.subtracted_dft_fields is not None:
+            self.sim.load_minus_flux_data(
+                self._monitor,
+                self.subtracted_dft_fields,
+            )
         return self._monitor
 
     def place_adjoint_source(self, dJ):
@@ -372,13 +391,15 @@ def __init__(
         Near2FarRegions: List[mp.Near2FarRegion],
         far_pts: List[mp.Vector3],
         decimation_factor: Optional[int] = 0,
+        norm_near_fields: Optional[NearToFarData] = None,
     ):
         """ """
         super().__init__(sim)
         self.Near2FarRegions = Near2FarRegions
         self.far_pts = far_pts  # list of far pts
         self._nfar_pts = len(far_pts)
         self.decimation_factor = decimation_factor
+        self.norm_near_fields = norm_near_fields
 
     def register_monitors(self, frequencies):
         self._frequencies = np.asarray(frequencies)
@@ -387,6 +408,11 @@ def register_monitors(self, frequencies):
             *self.Near2FarRegions,
             decimation_factor=self.decimation_factor,
         )
+        if self.norm_near_fields is not None:
+            self.sim.load_minus_near2far_data(
+                self._monitor,
+                self.norm_near_fields,
+            )
         return self._monitor
 
     def place_adjoint_source(self, dJ):

python/tests/test_adjoint_solver.py

@@ -137,6 +137,30 @@ def adjoint_solver(
 
         if mon_type.name == "EIGENMODE":
             if len(frequencies) == 1:
+                if mat2 is None:
+                    # Compute the incident fields of the mode source
+                    # in the straight waveguide for use as normalization
+                    # of the reflectance (S11) measurement.
+                    ref_sim = mp.Simulation(
+                        resolution=self.resolution,
+                        cell_size=self.cell_size,
+                        boundary_layers=self.pml_xy,
+                        sources=self.mode_source,
+                        geometry=self.waveguide_geometry,
+                    )
+                    dft_mon = ref_sim.add_mode_monitor(
+                        frequencies,
+                        mp.ModeRegion(
+                            center=mp.Vector3(-0.5 * self.sxy + self.dpml),
+                            size=mp.Vector3(0, self.sxy - 2 * self.dpml, 0),
+                        ),
+                        yee_grid=True,
+                    )
+                    ref_sim.run(until_after_sources=20)
+                    subtracted_dft_fields = ref_sim.get_flux_data(dft_mon)
+                else:
+                    subtracted_dft_fields = None
+
                 obj_list = [
                     mpa.EigenmodeCoefficient(
                         sim,
@@ -147,6 +171,7 @@ def adjoint_solver(
                         1,
                         forward=False,
                         eig_parity=self.eig_parity,
+                        subtracted_dft_fields=subtracted_dft_fields,
                     ),
                     mpa.EigenmodeCoefficient(
                         sim,