Revamped add fields::dft_fields_norm2

python/adjoint/optimization_problem.py

@@ -69,12 +69,10 @@ def __init__(
         fcen=None,
         df=None,
         nf=None,
-        decay_dt=50,
-        decay_fields=[mp.Ez],
-        decay_by=1e-6,
+        decay_by=1e-11,
         decimation_factor=0,
         minimum_run_time=0,
-        maximum_run_time=None,
+        maximum_run_time=None
     ):
 
         self.sim = simulation
@@ -130,8 +128,6 @@ def __init__(
                     2))  # index of center frequency
 
         self.decay_by = decay_by
-        self.decay_fields = decay_fields
-        self.decay_dt = decay_dt
         self.decimation_factor = decimation_factor
         self.minimum_run_time = minimum_run_time
         self.maximum_run_time = maximum_run_time
@@ -242,16 +238,10 @@ def forward_run(self):
         self.prepare_forward_run()
 
         # Forward run
-        self.sim.run(until_after_sources=stop_when_dft_decayed(
-            self.sim,
-            self.design_region_monitors,
-            self.decay_dt,
-            self.decay_fields,
-            self.fcen_idx,
+        self.sim.run(until_after_sources=mp.stop_when_dft_decayed(
             self.decay_by,
-            True,
             self.minimum_run_time,
-            self.maximum_run_time,
+            self.maximum_run_time
         ))
 
         # record objective quantities from user specified monitors
@@ -322,16 +312,10 @@ def adjoint_run(self):
             ]
 
             # Adjoint run
-            self.sim.run(until_after_sources=stop_when_dft_decayed(
-                self.sim,
-                self.design_region_monitors,
-                self.decay_dt,
-                self.decay_fields,
-                self.fcen_idx,
+            self.sim.run(until_after_sources=mp.stop_when_dft_decayed(
                 self.decay_by,
-                True,
                 self.minimum_run_time,
-                self.maximum_run_time,
+                self.maximum_run_time
             ))
 
             # Store adjoint fields for each design set of design variables in array (x,y,z,field_components,frequencies)
@@ -447,16 +431,10 @@ def calculate_fd_gradient(
                 self.forward_monitors.append(
                     m.register_monitors(self.frequencies))
 
-            self.sim.run(until_after_sources=stop_when_dft_decayed(
-                self.sim,
-                self.forward_monitors,
-                self.decay_dt,
-                self.decay_fields,
-                self.fcen_idx,
+            self.sim.run(until_after_sources=mp.stop_when_dft_decayed(
                 self.decay_by,
-                True,
                 self.minimum_run_time,
-                self.maximum_run_time,
+                self.maximum_run_time
             ))
 
             # record final objective function value
@@ -482,16 +460,10 @@ def calculate_fd_gradient(
                     m.register_monitors(self.frequencies))
 
             # add monitor used to track dft convergence
-            self.sim.run(until_after_sources=stop_when_dft_decayed(
-                self.sim,
-                self.forward_monitors,
-                self.decay_dt,
-                self.decay_fields,
-                self.fcen_idx,
+            self.sim.run(until_after_sources=mp.stop_when_dft_decayed(
                 self.decay_by,
-                True,
                 self.minimum_run_time,
-                self.maximum_run_time,
+                self.maximum_run_time
             ))
 
             # record final objective function value
@@ -548,93 +520,6 @@ def plot2D(self, init_opt=False, **kwargs):
 
         self.sim.plot2D(**kwargs)
 
-
-def stop_when_dft_decayed(
-    simob,
-    mon,
-    dt,
-    c,
-    fcen_idx,
-    decay_by,
-    yee_grid=False,
-    minimum_run_time=0,
-    maximum_run_time=None,
-):
-    '''Step function that monitors the relative change in DFT fields for a list of monitors.
-
-    mon ............. a list of monitors
-    c ............... a list of components to monitor
-
-    '''
-    # get monitor dft output array dimensions
-    dims = []
-    for m in mon:
-        ci_dims = []
-        for ci in c:
-            comp = ci if yee_grid else mp.Dielectric
-            ci_dims += [simob.get_array_slice_dimensions(comp, vol=m.where)[0]]
-        dims.append(ci_dims)
-
-    # Record data in closure so that we can persitently edit
-    closure = {'previous_fields': [[None for di in d] for d in dims], 't0': 0}
-
-    def _stop(sim):
-        if sim.round_time() <= dt + closure['t0']:
-            return False
-        elif maximum_run_time and sim.round_time() > maximum_run_time:
-            return True
-        else:
-            previous_fields = closure['previous_fields']
-
-            # Pull all the relevant frequency and spatial dft points
-            relative_change = []
-            current_fields = [[0 for di in d] for d in dims]
-            for mi, m in enumerate(mon):
-                for ic, cc in enumerate(c):
-                    if isinstance(m, mp.DftFlux):
-                        current_fields[mi][ic] = mp.get_fluxes(m)[fcen_idx]
-                    elif isinstance(m, mp.DftFields):
-                        current_fields[mi][ic] = atleast_3d(
-                            sim.get_dft_array(m, cc, fcen_idx))
-                    elif isinstance(m, mp.simulation.DftNear2Far):
-                        current_fields[mi][ic] = atleast_3d(
-                            sim.get_dft_array(m, cc, fcen_idx))
-                    else:
-                        raise TypeError(
-                            "Monitor of type {} not supported".format(type(m)))
-
-                    if previous_fields[mi][ic] is not None:
-                        if np.sum(
-                                np.abs(previous_fields[mi][ic] -
-                                       current_fields[mi][ic])) == 0:
-                            relative_change.append(0)
-                        elif np.all(np.abs(previous_fields[mi][ic])):
-                            relative_change_raw = np.abs(
-                                previous_fields[mi][ic] -
-                                current_fields[mi][ic]) / np.abs(
-                                    previous_fields[mi][ic])
-                            relative_change.append(
-                                np.mean(relative_change_raw.flatten())
-                            )  # average across space and frequency
-                        else:
-                            relative_change.append(1)
-                    else:
-                        relative_change.append(1)
-
-            relative_change = np.mean(
-                relative_change)  # average across monitors
-            closure['previous_fields'] = current_fields
-            closure['t0'] = sim.round_time()
-
-            if mp.verbosity > 0:
-                fmt = "DFT decay(t = {0:1.1f}): {1:0.4e}"
-                print(fmt.format(sim.meep_time(), np.real(relative_change)))
-            return relative_change <= decay_by and sim.round_time(
-            ) >= minimum_run_time
-
-    return _stop
-
-
 def atleast_3d(*arys):
     from numpy import array, asanyarray, newaxis
     '''

python/meep.i

@@ -1802,6 +1802,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
         scale_near2far_fields,
         stop_after_walltime,
         stop_on_interrupt,
+        stop_when_dft_decayed,
         stop_when_fields_decayed,
         synchronized_magnetic,
         to_appended,

python/simulation.py

@@ -4425,6 +4425,44 @@ def _stop(sim):
 
     return _stop
 
+def stop_when_dft_decayed(tol, minimum_run_time=0, maximum_run_time=None):
+    """
+    Return a `condition` function, suitable for passing to `sim.run()` as the `until`
+    or `until_after_sources` parameter, that checks all of the Simulation's dft objects
+    every `dt` timesteps, and stops the simulation once all the components and frequencies
+    of *every* dft object have decayed by at least some tolerance, `tol`. The routine 
+    calculates the optimal `dt` to use based on the frequency content in the DFT monitors.
+
+    Optionally the user can specify a minimum run time (`minimum_run_time`) or a maximum
+    run time (`maximum_run_time`).
+    """
+    # Record data in closure so that we can persitently edit
+    closure = {'previous_fields':0, 't0':0, 'dt':0, 'maxchange':0}
+    def _stop(_sim):
+        if _sim.fields.t == 0:
+            closure['dt'] = max(1/_sim.fields.dft_maxfreq()/_sim.fields.dt,_sim.fields.min_decimation())
+        if maximum_run_time and _sim.round_time() > maximum_run_time:
+            return True
+        elif _sim.fields.t <= closure['dt'] + closure['t0']:
+            return False
+        else:
+            previous_fields = closure['previous_fields']
+            current_fields  = _sim.fields.dft_norm()
+            change = np.abs(previous_fields-current_fields)
+            closure['maxchange'] = max(closure['maxchange'],change)
+
+            if previous_fields == 0:
+                closure['previous_fields'] = current_fields
+                return False
+
+            closure['previous_fields'] = current_fields
+            closure['t0'] = _sim.fields.t
+            if mp.verbosity > 0:
+                fmt = "DFT decay(t = {0:1.1f}): {1:0.4e}"
+                print(fmt.format(_sim.meep_time(), np.real(change/closure['maxchange'])))
+            return (change/closure['maxchange']) <= tol and _sim.round_time() >= minimum_run_time
+
+    return _stop
 
 def combine_step_funcs(*step_funcs):
     """

python/tests/test_adjoint_cyl.py

@@ -109,8 +109,6 @@ def J(alpha):
         fcen=fcen,
         df = 0,
         nf = 1,
-        decay_fields=[mp.Er],
-        decay_by=1e-4,
         maximum_run_time=1200)
 
     f, dJ_du = opt([design_params])

python/tests/test_adjoint_solver.py

@@ -164,7 +164,6 @@ def J(mode_mon):
         objective_arguments = obj_list,
         design_regions = [matgrid_region],
         frequencies=frequencies,
-        decay_fields=[mp.Ez],
         decimation_factor=10)
 
     f, dJ_du = opt([design_params])

src/dft.cpp

@@ -283,6 +283,81 @@ void dft_chunk::update_dft(double time) {
   }
 }
 
+/* Return the L2 norm of the DFTs themselves.  This is useful
+   to check whether the simulation is finished (whether all relevant fields have decayed).
+   (Collective operation.) */
+double fields::dft_norm() {
+  am_now_working_on(Other);
+  double sum = 0.0;
+  for (int i = 0; i < num_chunks; i++)
+    if (chunks[i]->is_mine()) sum += chunks[i]->dft_norm2();
+  finished_working();
+  return std::sqrt(sum_to_all(sum));
+}
+
+double fields_chunk::dft_norm2() const {
+  double sum = 0.0;
+  for (dft_chunk *cur = dft_chunks; cur; cur = cur->next_in_chunk)
+    sum += cur->norm2();
+  return sum;
+}
+
+static double sqr(std::complex<realnum> x) { return (x*std::conj(x)).real(); }
+
+double dft_chunk::norm2() const {
+  if (!fc->f[c][0]) return 0.0;
+  int numcmp = fc->f[c][1] ? 2 : 1;
+  double sum = 0.0;
+  size_t idx_dft = 0;
+  const int Nomega = omega.size();
+  LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
+    for (int i = 0; i < Nomega; ++i)
+        sum += sqr(dft[Nomega * idx_dft + i]);
+    idx_dft++;
+  }
+  return sum;
+}
+
+// return the minimum decimation factor across
+// all dft regions
+int fields::min_decimation() const {
+  int mindec = std::numeric_limits<int>::max();
+  for (int i = 0; i < num_chunks; i++)
+    if (chunks[i]->is_mine())
+      mindec = std::min(mindec, chunks[i]->min_decimation());
+  return min_to_all(mindec);
+}
+
+int fields_chunk::min_decimation() const {
+  int mindec = std::numeric_limits<int>::max();
+  for (dft_chunk *cur = dft_chunks; cur; cur = cur->next_in_chunk)
+    mindec = std::min(mindec, cur->get_decimation_factor());
+  return mindec;
+}
+
+// return the maximum abs(freq) over all DFT chunks
+double fields::dft_maxfreq() const {
+  double maxfreq = 0;
+  for (int i = 0; i < num_chunks; i++)
+    if (chunks[i]->is_mine())
+      maxfreq = std::max(maxfreq, chunks[i]->dft_maxfreq());
+  return max_to_all(maxfreq);
+}
+
+double fields_chunk::dft_maxfreq() const {
+  double maxomega = 0;
+  for (dft_chunk *cur = dft_chunks; cur; cur = cur->next_in_chunk)
+    maxomega = std::max(maxomega, cur->maxomega());
+  return maxomega / (2*meep::pi);
+}
+
+double dft_chunk::maxomega() const {
+  double maxomega = 0;
+  for (const auto& o : omega)
+    maxomega = std::max(maxomega, std::abs(o));
+  return maxomega;
+}
+
 void dft_chunk::scale_dft(complex<double> scale) {
   for (size_t i = 0; i < N * omega.size(); ++i)
     dft[i] *= scale;

src/meep.hpp

@@ -1114,6 +1114,8 @@ class dft_chunk {
   ~dft_chunk();
 
   void update_dft(double time);
+  double norm2() const;
+  double maxomega() const;
 
   void scale_dft(std::complex<double> scale);
 
@@ -1565,6 +1567,9 @@ class fields_chunk {
   void initialize_with_nth_tm(int n, double kz);
   // dft.cpp
   void update_dfts(double timeE, double timeH, int current_step);
+  double dft_norm2() const;
+  double dft_maxfreq() const;
+  int min_decimation() const;
 
   void changing_structure();
 };
@@ -1979,6 +1984,10 @@ class fields {
   dft_chunk *add_dft(const volume_list *where, const std::vector<double> &freq,
                      bool include_dV = true);
   void update_dfts();
+  double dft_norm();
+  double dft_maxfreq() const;
+  int min_decimation() const;
+
   dft_flux add_dft_flux(const volume_list *where, const double *freq, size_t Nfreq,
                         bool use_symmetry = true, bool centered_grid = true,
                         int decimation_factor = 0);

src/meep/mympi.hpp

@@ -72,6 +72,7 @@ bool broadcast(int from, bool);
 double max_to_master(double); // Only returns the correct value to proc 0.
 double max_to_all(double);
 int max_to_all(int);
+int min_to_all(int);
 float sum_to_master(float);   // Only returns the correct value to proc 0.
 double sum_to_master(double); // Only returns the correct value to proc 0.
 double sum_to_all(double);

src/mympi.cpp

@@ -387,6 +387,14 @@ int max_to_all(int in) {
   return out;
 }
 
+int min_to_all(int in) {
+  int out = in;
+#ifdef HAVE_MPI
+  MPI_Allreduce(&in, &out, 1, MPI_INT, MPI_MIN, mycomm);
+#endif
+  return out;
+}
+
 ivec max_to_all(const ivec &pt) {
   int in[5], out[5];
   for (int i = 0; i < 5; ++i)