enable automatic DFT decimation by default (#1732)

* automatic DFT decimation

* disable automatic decimation for certain tests

* fixes

* disable decimation for more tests

* compute the maximum frequency of the DFT monitor rather than use the last array element

* replace get_width with get_fwidth and use analytic formula for gaussian_src_time

* add missing 2pi factor to gaussian_src_time::get_fwidth

* add factor of 2 to gaussian_src_time::get_fwidth(double tol)

* various fixes

* disable automatic decimation for test_n2f_periodic.py

* switch ceil to floor when setting decimation_factor

doc/docs/Python_User_Interface.md

@@ -107,6 +107,7 @@ def __init__(self,
              progress_interval=4,
              subpixel_tol=0.0001,
              subpixel_maxeval=100000,
+             loop_tile_base=0,
              ensure_periodicity=True,
              num_chunks=0,
              Courant=0.5,
@@ -709,7 +710,7 @@ of the field at that point.
 
 ```python
 def add_dft_fields(self, *args, **kwargs):
-def add_dft_fields(cs, fcen, df, nfreq, freq, where=None, center=None, size=None, yee_grid=False, decimation_factor=1):
+def add_dft_fields(cs, fcen, df, nfreq, freq, where=None, center=None, size=None, yee_grid=False, decimation_factor=0):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -723,13 +724,14 @@ The volume can also be specified via the `center` and `size` arguments. The
 default routine interpolates the Fourier-transformed fields at the center of each
 voxel within the specified volume. Alternatively, the exact Fourier-transformed
 fields evaluated at each corresponding Yee grid point is available by setting
-`yee_grid` to `True`. To reduce the memory-bandwidth burden of
-accumulating DFT fields, an integer `decimation_factor` >= 1 can be
-specified. DFT field values are updated every `decimation_factor`
-timesteps. Use this feature with care, as the decimated timeseries may be
-corrupted by [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
-The choice of decimation factor should take into account the properties of all sources
-in the simulation as well as the frequency range of the DFT field monitor.
+`yee_grid` to `True`. To reduce the memory-bandwidth burden of accumulating
+DFT fields, an integer `decimation_factor` can be specified for updating the DFT
+fields at every `decimation_factor` timesteps. If `decimation_factor` is 0 (the default),
+this value is automatically determined from the
+[Nyquist rate](https://en.wikipedia.org/wiki/Nyquist_rate) of the bandwidth-limited
+sources and this DFT monitor. It can be turned off by setting it to 1. Use this feature
+with care, as the decimated timeseries may be corrupted by
+[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
 
 </div>
 
@@ -1091,7 +1093,7 @@ Given a bunch of [`FluxRegion`](#fluxregion) objects, you can tell Meep to accum
 
 ```python
 def add_flux(self, *args, **kwargs):
-def add_flux(fcen, df, nfreq, freq, FluxRegions, decimation_factor=1):
+def add_flux(fcen, df, nfreq, freq, FluxRegions, decimation_factor=0):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -1103,11 +1105,14 @@ field Fourier transforms for `nfreq` equally spaced frequencies covering the
 frequency range `fcen-df/2` to `fcen+df/2` or an array/list `freq` for arbitrarily
 spaced frequencies. Return a *flux object*, which you can pass to the functions
 below to get the flux spectrum, etcetera. To reduce the memory-bandwidth burden of
-accumulating DFT fields, an integer `decimation_factor` >= 1 can be
-specified. DFT field values are updated every `decimation_factor`
-timesteps. Use this feature with care, as the decimated timeseries may be
-corrupted by [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
-The choice of decimation factor should take into account the properties of all sources
+accumulating DFT fields, an integer `decimation_factor` can be specified for updating the DFT
+fields at every `decimation_factor` timesteps. If `decimation_factor` is 0 (the default),
+this value is automatically determined from the
+[Nyquist rate](https://en.wikipedia.org/wiki/Nyquist_rate) of the bandwidth-limited
+sources and this DFT monitor. It can be turned off by setting it to 1. Use this feature
+with care, as the decimated timeseries may be corrupted by
+[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. The choice
+of decimation factor should take into account the properties of all sources
 in the simulation as well as the frequency range of the DFT field monitor.
 
 </div>
@@ -1355,7 +1360,7 @@ Technically, MPB computes `ωₙ(k)` and then inverts it with Newton's method to
 
 ```python
 def add_mode_monitor(self, *args, **kwargs):
-def add_mode_monitor(fcen, df, nfreq, freq, ModeRegions, decimation_factor=1):
+def add_mode_monitor(fcen, df, nfreq, freq, ModeRegions, decimation_factor=0):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -1485,7 +1490,7 @@ The usage is similar to the flux spectra: you define a set of [`EnergyRegion`](#
 
 ```python
 def add_energy(self, *args, **kwargs):
-def add_energy(fcen, df, nfreq, freq, EnergyRegions, decimation_factor=1):
+def add_energy(fcen, df, nfreq, freq, EnergyRegions, decimation_factor=0):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -1497,12 +1502,13 @@ field Fourier transforms for `nfreq` equally spaced frequencies covering the
 frequency range `fcen-df/2` to `fcen+df/2` or an array/list `freq` for arbitrarily
 spaced frequencies. Return an *energy object*, which you can pass to the functions
 below to get the energy spectrum, etcetera. To reduce the memory-bandwidth burden of
-accumulating DFT fields, an integer `decimation_factor` >= 1 can be
-specified. DFT field values are updated every `decimation_factor`
-timesteps. Use this feature with care, as the decimated timeseries may be
-corrupted by [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
-The choice of decimation factor should take into account the properties of all sources
-in the simulation as well as the frequency range of the DFT field monitor.
+accumulating DFT fields, an integer `decimation_factor` can be specified for updating the DFT
+fields at every `decimation_factor` timesteps. If `decimation_factor` is 0 (the default),
+this value is automatically determined from the
+[Nyquist rate](https://en.wikipedia.org/wiki/Nyquist_rate) of the bandwidth-limited
+sources and this DFT monitor. It can be turned off by setting it to 1. Use this feature
+with care, as the decimated timeseries may be corrupted by
+[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
 
 </div>
 
@@ -1717,7 +1723,7 @@ The usage is similar to the [flux spectra](Python_Tutorials/Basics.md#transmitta
 
 ```python
 def add_force(self, *args, **kwargs):
-def add_force(fcen, df, nfreq, freq, ForceRegions, decimation_factor=1):
+def add_force(fcen, df, nfreq, freq, ForceRegions, decimation_factor=0):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -1729,12 +1735,13 @@ field Fourier transforms for `nfreq` equally spaced frequencies covering the
 frequency range `fcen-df/2` to `fcen+df/2` or an array/list `freq` for arbitrarily
 spaced frequencies. Return a `force`object, which you can pass to the functions
 below to get the force spectrum, etcetera. To reduce the memory-bandwidth burden of
-accumulating DFT fields, an integer `decimation_factor` >= 1 can be
-specified. DFT field values are updated every `decimation_factor`
-timesteps. Use this feature with care, as the decimated timeseries may be
-corrupted by [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
-The choice of decimation factor should take into account the properties of all sources
-in the simulation as well as the frequency range of the DFT field monitor.
+accumulating DFT fields, an integer `decimation_factor` can be specified for updating the DFT
+fields at every `decimation_factor` timesteps. If `decimation_factor` is 0 (the default),
+this value is automatically determined from the
+[Nyquist rate](https://en.wikipedia.org/wiki/Nyquist_rate) of the bandwidth-limited
+sources and this DFT monitor. It can be turned off by setting it to 1. Use this feature
+with care, as the decimated timeseries may be corrupted by
+[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
 
 </div>
 
@@ -2002,7 +2009,7 @@ There are three steps to using the near-to-far-field feature: first, define the
 
 ```python
 def add_near2far(self, *args, **kwargs):
-def add_near2far(fcen, df, nfreq, freq, Near2FarRegions, nperiods=1, decimation_factor=1):
+def add_near2far(fcen, df, nfreq, freq, Near2FarRegions, nperiods=1, decimation_factor=0):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -2014,12 +2021,13 @@ appropriate field Fourier transforms for `nfreq` equally spaced frequencies
 covering the frequency range `fcen-df/2` to `fcen+df/2` or an array/list `freq`
 for arbitrarily spaced frequencies. Return a `near2far` object, which you can pass
 to the functions below to get the far fields. To reduce the memory-bandwidth burden of
-accumulating DFT fields, an integer `decimation_factor` >= 1 can be
-specified. DFT field values are updated every `decimation_factor`
-timesteps. Use this feature with care, as the decimated timeseries may be
-corrupted by [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
-The choice of decimation factor should take into account the properties of all sources
-in the simulation as well as the frequency range of the DFT field monitor.
+accumulating DFT fields, an integer `decimation_factor` can be specified for updating the DFT
+fields at every `decimation_factor` timesteps. If `decimation_factor` is 0 (the default),
+this value is automatically determined from the
+[Nyquist rate](https://en.wikipedia.org/wiki/Nyquist_rate) of the bandwidth-limited
+sources and this DFT monitor. It can be turned off by setting it to 1. Use this feature
+with care, as the decimated timeseries may be corrupted by
+[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
 
 </div>
 
@@ -4262,8 +4270,7 @@ $\beta$ where $\beta=+\infty$ gives an unsmoothed, discontinuous interface. The
 +\tanh(\beta\times(u-\eta)))/(\tanh(\beta\times\eta)+\tanh(\beta\times(1-\eta)))$ involving the parameters `beta`
 ($\beta$: bias or "smoothness" of the turn on) and `eta` ($\eta$: offset for erosion/dilation). The level set provides a general
 approach for defining a *discontinuous* function from otherwise continuously varying (via the bilinear interpolation) grid values.
-The subpixel smoothing is based on an adaptive quadrature scheme with properties `subpixel_maxeval` and `subpixel_tol` which
-can be specified using the `Simulation` constructor.
+Subpixel smoothing is fast and accurate because it exploits an analytic formulation for level-set functions.
 
 It is possible to overlap any number of different `MaterialGrid`s. This can be useful for defining grids which are symmetric
 (e.g., mirror, rotation). One way to set this up is by overlapping a given `MaterialGrid` object with a symmetrized copy of

python/adjoint/objective.py

@@ -121,7 +121,7 @@ def __init__(self,
                  mode,
                  forward=True,
                  kpoint_func=None,
-                 decimation_factor=1,
+                 decimation_factor=0,
                  **kwargs):
         super().__init__(sim)
         self.volume = volume
@@ -206,7 +206,7 @@ def __call__(self):
 
 
 class FourierFields(ObjectiveQuantity):
-    def __init__(self, sim, volume, component, decimation_factor=1):
+    def __init__(self, sim, volume, component, decimation_factor=0):
         super().__init__(sim)
         self.volume = volume
         self.component = component
@@ -306,7 +306,7 @@ def __call__(self):
 
 
 class Near2FarFields(ObjectiveQuantity):
-    def __init__(self, sim, Near2FarRegions, far_pts, decimation_factor=1):
+    def __init__(self, sim, Near2FarRegions, far_pts, decimation_factor=0):
         super().__init__(sim)
         self.Near2FarRegions = Near2FarRegions
         self.far_pts = far_pts  #list of far pts

python/adjoint/optimization_problem.py

@@ -72,7 +72,7 @@ def __init__(
         decay_dt=50,
         decay_fields=[mp.Ez],
         decay_by=1e-6,
-        decimation_factor=1,
+        decimation_factor=0,
         minimum_run_time=0,
         maximum_run_time=None,
     ):

python/adjoint/utils.py

@@ -69,6 +69,7 @@ def install_design_region_monitors(
             frequencies,
             where=design_region.volume,
             yee_grid=True,
+            decimation_factor=0
         ) for design_region in design_regions
     ]
     return design_region_monitors

python/geom.py

@@ -547,7 +547,7 @@ def __init__(self,grid_size,medium1,medium2,weights=None,grid_type="U_DEFAULT",d
         +\\tanh(\\beta\\times(u-\\eta)))/(\\tanh(\\beta\\times\\eta)+\\tanh(\\beta\\times(1-\\eta)))$ involving the parameters `beta`
         ($\\beta$: bias or "smoothness" of the turn on) and `eta` ($\\eta$: offset for erosion/dilation). The level set provides a general
         approach for defining a *discontinuous* function from otherwise continuously varying (via the bilinear interpolation) grid values.
-        Subpixel smoothing is fast and accurate because it exploits an analytic framework for level-set functions.
+        Subpixel smoothing is fast and accurate because it exploits an analytic formulation for level-set functions.
 
         It is possible to overlap any number of different `MaterialGrid`s. This can be useful for defining grids which are symmetric
         (e.g., mirror, rotation). One way to set this up is by overlapping a given `MaterialGrid` object with a symmetrized copy of