diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index c539ea90b..33f59f655 100644 --- a/doc/docs/Python_User_Interface.md +++ b/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, @@ -715,7 +716,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): ```
@@ -729,13 +730,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.
@@ -1097,7 +1099,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): ```
@@ -1109,11 +1111,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.
@@ -1361,7 +1366,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): ```
@@ -1491,7 +1496,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): ```
@@ -1503,12 +1508,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.
@@ -1723,7 +1729,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): ```
@@ -1735,12 +1741,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.
@@ -2008,7 +2015,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): ```
@@ -2020,12 +2027,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.
@@ -4266,8 +4274,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 diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 802b5571a..e1d0159fb 100644 --- a/python/adjoint/objective.py +++ b/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 diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index 244103116..bb938cf45 100644 --- a/python/adjoint/optimization_problem.py +++ b/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, ): diff --git a/python/adjoint/utils.py b/python/adjoint/utils.py index 110a4c4c2..decbe9e11 100644 --- a/python/adjoint/utils.py +++ b/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 diff --git a/python/geom.py b/python/geom.py index 861b62d93..b061ad81e 100755 --- a/python/geom.py +++ b/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 diff --git a/python/simulation.py b/python/simulation.py index 62a8b54a4..34231ab6e 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -2413,7 +2413,7 @@ def _evaluate_dft_objects(self): def add_dft_fields(self, *args, **kwargs): """ - `add_dft_fields(cs, fcen, df, nfreq, freq, where=None, center=None, size=None, yee_grid=False, decimation_factor=1)` ##sig + `add_dft_fields(cs, fcen, df, nfreq, freq, where=None, center=None, size=None, yee_grid=False, decimation_factor=0)` ##sig Given a list of field components `cs`, compute the Fourier transform of these fields for `nfreq` equally spaced frequencies covering the frequency range @@ -2423,13 +2423,14 @@ def add_dft_fields(self, *args, **kwargs): 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. """ components = args[0] args = fix_dft_args(args, 1) @@ -2438,7 +2439,7 @@ def add_dft_fields(self, *args, **kwargs): center = kwargs.get('center', None) size = kwargs.get('size', None) yee_grid = kwargs.get('yee_grid', False) - decimation_factor = kwargs.get('decimation_factor', 1) + decimation_factor = kwargs.get('decimation_factor', 0) center_v3 = Vector3(*center) if center is not None else None size_v3 = Vector3(*size) if size is not None else None use_centered_grid = not yee_grid @@ -2487,7 +2488,7 @@ def get_dft_data(self, dft_chunk): def add_near2far(self, *args, **kwargs): """ - `add_near2far(fcen, df, nfreq, freq, Near2FarRegions, nperiods=1, decimation_factor=1)` ##sig + `add_near2far(fcen, df, nfreq, freq, Near2FarRegions, nperiods=1, decimation_factor=0)` ##sig Add a bunch of `Near2FarRegion`s to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the @@ -2495,18 +2496,19 @@ def add_near2far(self, *args, **kwargs): 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. """ args = fix_dft_args(args, 0) freq = args[0] near2fars = args[1:] nperiods = kwargs.get('nperiods', 1) - decimation_factor = kwargs.get('decimation_factor', 1) + decimation_factor = kwargs.get('decimation_factor', 0) n2f = DftNear2Far(self._add_near2far, [freq, nperiods, near2fars, decimation_factor]) self.dft_objects.append(n2f) return n2f @@ -2519,7 +2521,7 @@ def _add_near2far(self, freq, nperiods, near2fars, decimation_factor): def add_energy(self, *args, **kwargs): """ - `add_energy(fcen, df, nfreq, freq, EnergyRegions, decimation_factor=1)` ##sig + `add_energy(fcen, df, nfreq, freq, EnergyRegions, decimation_factor=0)` ##sig Add a bunch of `EnergyRegion`s to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the appropriate @@ -2527,17 +2529,18 @@ def add_energy(self, *args, **kwargs): 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. """ args = fix_dft_args(args, 0) freq = args[0] energys = args[1:] - decimation_factor = kwargs.get('decimation_factor', 1) + decimation_factor = kwargs.get('decimation_factor', 0) en = DftEnergy(self._add_energy, [freq, energys, decimation_factor]) self.dft_objects.append(en) return en @@ -2740,7 +2743,7 @@ def load_minus_near2far_data(self, near2far, n2fdata): def add_force(self, *args, **kwargs): """ - `add_force(fcen, df, nfreq, freq, ForceRegions, decimation_factor=1)` ##sig + `add_force(fcen, df, nfreq, freq, ForceRegions, decimation_factor=0)` ##sig Add a bunch of `ForceRegion`s to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the appropriate @@ -2748,17 +2751,18 @@ def add_force(self, *args, **kwargs): 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. """ args = fix_dft_args(args, 0) freq = args[0] forces = args[1:] - decimation_factor = kwargs.get('decimation_factor', 1) + decimation_factor = kwargs.get('decimation_factor', 0) force = DftForce(self._add_force, [freq, forces, decimation_factor]) self.dft_objects.append(force) return force @@ -2845,7 +2849,7 @@ def load_minus_force_data(self, force, fdata): def add_flux(self, *args, **kwargs): """ - `add_flux(fcen, df, nfreq, freq, FluxRegions, decimation_factor=1)` ##sig + `add_flux(fcen, df, nfreq, freq, FluxRegions, decimation_factor=0)` ##sig Add a bunch of `FluxRegion`s to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the appropriate @@ -2853,17 +2857,20 @@ def add_flux(self, *args, **kwargs): 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. """ args = fix_dft_args(args, 0) freq = args[0] fluxes = args[1:] - decimation_factor = kwargs.get('decimation_factor', 1) + decimation_factor = kwargs.get('decimation_factor', 0) flux = DftFlux(self._add_flux, [freq, fluxes, decimation_factor]) self.dft_objects.append(flux) return flux @@ -2876,14 +2883,14 @@ def _add_flux(self, freq, fluxes, decimation_factor): def add_mode_monitor(self, *args, **kwargs): """ - `add_mode_monitor(fcen, df, nfreq, freq, ModeRegions, decimation_factor=1)` ##sig + `add_mode_monitor(fcen, df, nfreq, freq, ModeRegions, decimation_factor=0)` ##sig Similar to `add_flux`, but for use with `get_eigenmode_coefficients`. """ args = fix_dft_args(args, 0) freq = args[0] fluxes = args[1:] - decimation_factor = kwargs.get('decimation_factor', 1) + decimation_factor = kwargs.get('decimation_factor', 0) yee_grid = kwargs.get("yee_grid", False) flux = DftFlux(self._add_mode_monitor, [freq, fluxes, yee_grid, decimation_factor]) self.dft_objects.append(flux) diff --git a/python/tests/test_3rd_harm_1d.py b/python/tests/test_3rd_harm_1d.py index f0bf9269f..87e9ee4e7 100644 --- a/python/tests/test_3rd_harm_1d.py +++ b/python/tests/test_3rd_harm_1d.py @@ -34,9 +34,9 @@ def setUp(self): dimensions=dimensions) fr = mp.FluxRegion(mp.Vector3(0, 0, (0.5 * self.sz) - self.dpml - 0.5)) - self.trans = self.sim.add_flux(0.5 * (fmin + fmax), fmax - fmin, nfreq, fr) - self.trans1 = self.sim.add_flux(fcen, 0, 1, fr) - self.trans3 = self.sim.add_flux(3 * fcen, 0, 1, fr) + self.trans = self.sim.add_flux(0.5 * (fmin + fmax), fmax - fmin, nfreq, fr, decimation_factor=1) + self.trans1 = self.sim.add_flux(fcen, 0, 1, fr, decimation_factor=1) + self.trans3 = self.sim.add_flux(3 * fcen, 0, 1, fr, decimation_factor=1) def test_3rd_harm_1d(self): diff --git a/python/tests/test_adjoint_jax.py b/python/tests/test_adjoint_jax.py index 346946b13..0cc927b3f 100644 --- a/python/tests/test_adjoint_jax.py +++ b/python/tests/test_adjoint_jax.py @@ -117,7 +117,8 @@ def build_straight_wg_simulation( simulation, mp.Volume(center=center, size=monitor_size), mode=1, - forward=True) for center in monitor_centers + forward=True, + decimation_factor=1) for center in monitor_centers ] return simulation, sources, monitors, design_regions, frequencies diff --git a/python/tests/test_bend_flux.py b/python/tests/test_bend_flux.py index 17dff7950..a21d7c71f 100644 --- a/python/tests/test_bend_flux.py +++ b/python/tests/test_bend_flux.py @@ -62,12 +62,13 @@ def init(self, no_bend=False, gdsii=False): else: fr = mp.FluxRegion(center=mp.Vector3(wvg_xcen, (sy / 2) - 1.5), size=mp.Vector3(w * 2, 0)) - self.trans = self.sim.add_flux(fcen, df, nfreq, fr) + self.trans = self.sim.add_flux(fcen, df, nfreq, fr, decimation_factor=1) self.trans_decimated = self.sim.add_flux(fcen, df, nfreq, fr, decimation_factor=5) refl_fr = mp.FluxRegion(center=mp.Vector3((-0.5 * sx) + 1.5, wvg_ycen), size=mp.Vector3(0, w * 2)) - self.refl = self.sim.add_flux(np.linspace(fcen-0.5*df,fcen+0.5*df,nfreq), refl_fr) + self.refl = self.sim.add_flux(np.linspace(fcen-0.5*df,fcen+0.5*df,nfreq), refl_fr, + decimation_factor=1) self.refl_decimated = self.sim.add_flux(np.linspace(fcen-0.5*df,fcen+0.5*df,nfreq), refl_fr, decimation_factor=10) diff --git a/python/tests/test_cavity_farfield.py b/python/tests/test_cavity_farfield.py index 861831c94..93de61f01 100644 --- a/python/tests/test_cavity_farfield.py +++ b/python/tests/test_cavity_farfield.py @@ -56,7 +56,8 @@ def run_test(self, nfreqs): size=mp.Vector3(0, d1), weight=-1.0), mp.Near2FarRegion(mp.Vector3(0.5 * sx - dpml, 0.5 * w + 0.5 * d1), - size=mp.Vector3(0, d1)) + size=mp.Vector3(0, d1)), + decimation_factor=1 ) sim.run(until=200) d2 = 20 diff --git a/python/tests/test_chunks.py b/python/tests/test_chunks.py index 886f57d4e..2af803b07 100644 --- a/python/tests/test_chunks.py +++ b/python/tests/test_chunks.py @@ -38,7 +38,7 @@ def test_chunks(self): rgt = mp.FluxRegion(center=mp.Vector3(+0.5*sxy-dpml,0), size=mp.Vector3(0,sxy-2*dpml), weight=+1.0) lft = mp.FluxRegion(center=mp.Vector3(-0.5*sxy+dpml,0), size=mp.Vector3(0,sxy-2*dpml), weight=-1.0) - tot_flux = sim.add_flux(fcen, 0, 1, top, bot, rgt, lft) + tot_flux = sim.add_flux(fcen, 0, 1, top, bot, rgt, lft, decimation_factor=1) sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(), 1e-5)) @@ -61,7 +61,7 @@ def test_chunks(self): sim.use_output_directory(self.temp_dir) - tot_flux = sim.add_flux(fcen, 0, 1, top, bot, rgt, lft) + tot_flux = sim.add_flux(fcen, 0, 1, top, bot, rgt, lft, decimation_factor=1) sim.load_minus_flux('tot_flux', tot_flux) diff --git a/python/tests/test_dft_energy.py b/python/tests/test_dft_energy.py index 856697a81..db43946ea 100644 --- a/python/tests/test_dft_energy.py +++ b/python/tests/test_dft_energy.py @@ -22,8 +22,10 @@ def test_dft_energy(self): sim = mp.Simulation(resolution=resolution, cell_size=cell, geometry=geom, boundary_layers=pml, sources=sources, symmetries=[mp.Mirror(direction=mp.Y)]) - flux = sim.add_flux(fsrc, 0, 1, mp.FluxRegion(center=mp.Vector3(3), size=mp.Vector3(y=5))) - energy = sim.add_energy(fsrc, 0, 1, mp.EnergyRegion(center=mp.Vector3(3), size=mp.Vector3(y=5))) + flux = sim.add_flux(fsrc, 0, 1, mp.FluxRegion(center=mp.Vector3(3), size=mp.Vector3(y=5)), + decimation_factor=1) + energy = sim.add_energy(fsrc, 0, 1, mp.EnergyRegion(center=mp.Vector3(3), size=mp.Vector3(y=5)), + decimation_factor=1) energy_decimated = sim.add_energy(fsrc, 0, 1, mp.EnergyRegion(center=mp.Vector3(3), size=mp.Vector3(y=5)), decimation_factor=10) diff --git a/python/tests/test_dft_fields.py b/python/tests/test_dft_fields.py index a66954abb..3b190adbe 100644 --- a/python/tests/test_dft_fields.py +++ b/python/tests/test_dft_fields.py @@ -102,7 +102,8 @@ def test_get_dft_array(self): def test_decimated_dft_fields_are_almost_equal_to_undecimated_fields(self): sim = self.init() sim.init_sim() - undecimated_field = sim.add_dft_fields([mp.Ez], self.fcen, 0, 1) + undecimated_field = sim.add_dft_fields([mp.Ez], self.fcen, 0, 1, + decimation_factor=1) decimated_field = sim.add_dft_fields([mp.Ez], self.fcen, 0, diff --git a/python/tests/test_divide_mpi_processes.py b/python/tests/test_divide_mpi_processes.py index ba7addf3c..574f0e834 100644 --- a/python/tests/test_divide_mpi_processes.py +++ b/python/tests/test_divide_mpi_processes.py @@ -1,5 +1,3 @@ -from __future__ import division - import unittest import meep as mp @@ -36,7 +34,8 @@ def test_divide_parallel_processes(self): mp.FluxRegion(mp.Vector3(y=0.5*sxy), size=mp.Vector3(sxy)), mp.FluxRegion(mp.Vector3(y=-0.5*sxy), size=mp.Vector3(sxy), weight=-1), mp.FluxRegion(mp.Vector3(0.5*sxy), size=mp.Vector3(y=sxy)), - mp.FluxRegion(mp.Vector3(-0.5*sxy), size=mp.Vector3(y=sxy), weight=-1)) + mp.FluxRegion(mp.Vector3(-0.5*sxy), size=mp.Vector3(y=sxy), weight=-1), + decimation_factor=1) self.sim.run(until_after_sources=30) diff --git a/python/tests/test_force.py b/python/tests/test_force.py index 136202d0d..a2588de45 100644 --- a/python/tests/test_force.py +++ b/python/tests/test_force.py @@ -20,7 +20,7 @@ def setUp(self): sources=[sources]) fr = mp.ForceRegion(mp.Vector3(y=1.27), direction=mp.Y, size=mp.Vector3(4.38)) - self.myforce = self.sim.add_force(fcen, 0, 1, fr) + self.myforce = self.sim.add_force(fcen, 0, 1, fr, decimation_factor=1) self.myforce_decimated = self.sim.add_force(fcen, 0, 1, fr, decimation_factor=10) def test_force(self): diff --git a/python/tests/test_geom.py b/python/tests/test_geom.py index a0fe5325e..4b2f3a974 100644 --- a/python/tests/test_geom.py +++ b/python/tests/test_geom.py @@ -348,12 +348,12 @@ def check_warnings(sim, should_warn=True): sim = mp.Simulation(cell_size=cell_size, resolution=resolution, sources=valid_sources[0], extra_materials=[mat]) fregion = mp.FluxRegion(center=mp.Vector3(0, 1), size=mp.Vector3(2, 2), direction=mp.X) - sim.add_flux(18, 6, 2, fregion) + sim.add_flux(18, 6, 2, fregion, decimation_factor=1) check_warnings(sim) # Invalid geometry material sim = mp.Simulation(cell_size=cell_size, resolution=resolution, sources=valid_sources[0], geometry=geom) - sim.add_flux(18, 6, 2, fregion) + sim.add_flux(18, 6, 2, fregion, decimation_factor=1) check_warnings(sim) def test_transform(self): diff --git a/python/tests/test_holey_wvg_cavity.py b/python/tests/test_holey_wvg_cavity.py index 5b486273f..a66998d28 100644 --- a/python/tests/test_holey_wvg_cavity.py +++ b/python/tests/test_holey_wvg_cavity.py @@ -134,7 +134,8 @@ def test_transmission_spectrum(self): freg = mp.FluxRegion(center=mp.Vector3((0.5 * self.sx) - self.dpml - 0.5), size=mp.Vector3(0, 2 * self.w)) - trans = self.sim.add_flux(self.fcen, self.df, self.nfreq, freg) + trans = self.sim.add_flux(self.fcen, self.df, self.nfreq, freg, + decimation_factor=1) self.sim.run( until_after_sources=mp.stop_when_fields_decayed( diff --git a/python/tests/test_mode_decomposition.py b/python/tests/test_mode_decomposition.py index b6abcd0e4..6c400cf56 100644 --- a/python/tests/test_mode_decomposition.py +++ b/python/tests/test_mode_decomposition.py @@ -116,7 +116,8 @@ def test_oblique_waveguide_backward_mode(self): mode = sim.add_mode_monitor(fcen, 0, 1, mp.FluxRegion(center=mp.Vector3(-0.5*sxy+dpml,0,0), - size=mp.Vector3(0,sxy,0))) + size=mp.Vector3(0,sxy,0)), + decimation_factor=1) mode_decimated = sim.add_mode_monitor(fcen, 0, 1, mp.FluxRegion(center=mp.Vector3(-0.5*sxy+dpml,0,0), size=mp.Vector3(0,sxy,0)), diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py index 89acd03a4..209dd16eb 100644 --- a/python/tests/test_refl_angular.py +++ b/python/tests/test_refl_angular.py @@ -41,7 +41,7 @@ def test_refl_angular(self): resolution=resolution) refl_fr = mp.FluxRegion(center=mp.Vector3(z=-0.25 * sz)) - refl = sim.add_flux(fcen, df, nfreq, refl_fr) + refl = sim.add_flux(fcen, df, nfreq, refl_fr, decimation_factor=1) sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(z=-0.5 * sz + dpml), 1e-9)) @@ -59,7 +59,7 @@ def test_refl_angular(self): dimensions=dimensions, resolution=resolution) - refl = sim.add_flux(fcen, df, nfreq, refl_fr) + refl = sim.add_flux(fcen, df, nfreq, refl_fr, decimation_factor=1) sim.load_minus_flux_data(refl, empty_data) sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(z=-0.5 * sz + dpml), 1e-9)) diff --git a/python/tests/test_visualization.py b/python/tests/test_visualization.py index 4400955e4..e2b9d2846 100644 --- a/python/tests/test_visualization.py +++ b/python/tests/test_visualization.py @@ -82,10 +82,10 @@ def setup_sim(zDim=0): resolution = 10 sim = mp.Simulation(cell_size=cell, - boundary_layers=pml_layers, - geometry=geometry, - sources=sources, - resolution=resolution) + boundary_layers=pml_layers, + geometry=geometry, + sources=sources, + resolution=resolution) # Line monitor sim.add_flux(1,0,1,mp.FluxRegion(center=mp.Vector3(5,0,0),size=mp.Vector3(0,4,4), direction=mp.X)) diff --git a/src/dft.cpp b/src/dft.cpp index cf5177caa..0ecf4825f 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -174,7 +174,26 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, dft_chunk_data data; data.c = c; data.vc = vc; + + if (decimation_factor == 0) { + double tol = 1e-7; + double src_freq_max = 0; + for (src_time *s = sources; s; s = s->next) { + if (s->get_fwidth(tol) == 0) + decimation_factor = 1; + else + src_freq_max = std::max(src_freq_max, std::abs(s->frequency().real())+0.5*s->get_fwidth(tol)); + } + double freq_max = 0; + for (size_t i = 0; i < Nfreq; ++i) + freq_max = std::max(freq_max, std::abs(freq[i])); + if ((freq_max > 0) && (src_freq_max > 0)) + decimation_factor = std::max(1, int(std::floor(1/(dt*(freq_max + src_freq_max))))); + else + decimation_factor = 1; + } data.decimation_factor = decimation_factor; + data.omega.resize(Nfreq); for (size_t i = 0; i < Nfreq; ++i) data.omega[i] = 2 * pi * freq[i]; diff --git a/src/meep.hpp b/src/meep.hpp index aeb28211e..7dc9fdf93 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -973,6 +973,7 @@ class src_time { return 1; } virtual std::complex frequency() const { return 0.0; } + virtual double get_fwidth(double tol) const { return 0.0; } virtual void set_frequency(std::complex f) { (void)f; } private: @@ -994,6 +995,7 @@ class gaussian_src_time : public src_time { virtual src_time *clone() const { return new gaussian_src_time(*this); } virtual bool is_equal(const src_time &t) const; virtual std::complex frequency() const { return freq; } + virtual double get_fwidth(double tol) const; virtual void set_frequency(std::complex f) { freq = real(f); } std::complex fourier_transform(const double f); @@ -1014,6 +1016,7 @@ class continuous_src_time : public src_time { virtual src_time *clone() const { return new continuous_src_time(*this); } virtual bool is_equal(const src_time &t) const; virtual std::complex frequency() const { return freq; } + virtual double get_fwidth(double tol) const { return 0.0; }; virtual void set_frequency(std::complex f) { freq = f; } private: @@ -1046,6 +1049,7 @@ class custom_src_time : public src_time { virtual src_time *clone() const { return new custom_src_time(*this); } virtual bool is_equal(const src_time &t) const; virtual std::complex frequency() const { return freq; } + virtual double get_fwidth(double tol) const { return 0.0; }; virtual void set_frequency(std::complex f) { freq = f; } private: @@ -1925,7 +1929,7 @@ class fields { std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, - int vc = 0, int decimation_factor = 1) { + int vc = 0, int decimation_factor = 0) { return add_dft(c, where, linspace(freq_min, freq_max, Nfreq), include_dV_and_interp_weights, stored_weight, chunk_next, sqrt_dV_and_interp_weights, extra_weight, use_centered_grid, vc, decimation_factor); @@ -1935,13 +1939,13 @@ class fields { std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, - int vc = 0, int decimation_factor = 1); + int vc = 0, int decimation_factor = 0); dft_chunk *add_dft(component c, const volume &where, const std::vector& freq, bool include_dV_and_interp_weights = true, std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, - int vc = 0, int decimation_factor = 1) { + int vc = 0, int decimation_factor = 0) { return add_dft(c, where, freq.data(), freq.size(), include_dV_and_interp_weights, stored_weight, chunk_next, sqrt_dV_and_interp_weights, extra_weight, use_centered_grid, vc, decimation_factor); @@ -1962,34 +1966,34 @@ class fields { void update_dfts(); 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 = 1); + int decimation_factor = 0); dft_flux add_dft_flux(const volume_list *where, const std::vector &freq, - int decimation_factor = 1, bool use_symmetry = true, + int decimation_factor = 0, bool use_symmetry = true, bool centered_grid = true) { return add_dft_flux(where, freq.data(), freq.size(), use_symmetry, centered_grid, decimation_factor); } dft_flux add_dft_flux(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true, bool centered_grid = true, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_flux(where, linspace(freq_min, freq_max, Nfreq), use_symmetry, centered_grid, decimation_factor); } dft_flux add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true, bool centered_grid = true, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_flux(d, where, linspace(freq_min, freq_max, Nfreq), use_symmetry, centered_grid, decimation_factor); } dft_flux add_dft_flux(direction d, const volume &where, const std::vector &freq, bool use_symmetry = true, bool centered_grid = true, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_flux(d, where, freq.data(), freq.size(), use_symmetry, centered_grid, decimation_factor); } dft_flux add_dft_flux(direction d, const volume &where, const double *freq, size_t Nfreq, bool use_symmetry = true, bool centered_grid = true, - int decimation_factor = 1); + int decimation_factor = 0); dft_flux add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq); dft_flux add_dft_flux_box(const volume &where, const std::vector &freq); dft_flux add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq); @@ -1997,33 +2001,33 @@ class fields { // a "mode monitor" is just a dft_flux with symmetry reduction turned off. dft_flux add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, - int Nfreq, bool centered_grid = true, int decimation_factor = 1) { + int Nfreq, bool centered_grid = true, int decimation_factor = 0) { return add_mode_monitor(d, where, linspace(freq_min, freq_max, Nfreq), centered_grid, decimation_factor); } dft_flux add_mode_monitor(direction d, const volume &where, const std::vector &freq, - bool centered_grid = true, int decimation_factor = 1) { + bool centered_grid = true, int decimation_factor = 0) { return add_mode_monitor(d, where, freq.data(), freq.size(), centered_grid, decimation_factor); } dft_flux add_mode_monitor(direction d, const volume &where, const double *freq, size_t Nfreq, - bool centered_grid = true, int decimation_factor = 1); + bool centered_grid = true, int decimation_factor = 0); dft_fields add_dft_fields(component *components, int num_components, const volume where, double freq_min, double freq_max, int Nfreq, - bool use_centered_grid = true, int decimation_factor = 1) { + bool use_centered_grid = true, int decimation_factor = 0) { return add_dft_fields(components, num_components, where, linspace(freq_min, freq_max, Nfreq), use_centered_grid, decimation_factor); } dft_fields add_dft_fields(component *components, int num_components, const volume where, const std::vector &freq, bool use_centered_grid = true, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_fields(components, num_components, where, freq.data(), freq.size(), use_centered_grid, decimation_factor); } dft_fields add_dft_fields(component *components, int num_components, const volume where, const double *freq, size_t Nfreq, bool use_centered_grid = true, - int decimation_factor = 1); + int decimation_factor = 0); /********************************************************/ /* process_dft_component is an intermediate-level */ @@ -2069,41 +2073,41 @@ class fields { std::complex overlaps[2]); dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_energy(where, linspace(freq_min, freq_max, Nfreq), decimation_factor); } dft_energy add_dft_energy(const volume_list *where, const std::vector &freq, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_energy(where, freq.data(), freq.size(), decimation_factor); } dft_energy add_dft_energy(const volume_list *where, const double *freq, size_t Nfreq, - int decimation_factor = 1); + int decimation_factor = 0); // stress.cpp dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_force(where, linspace(freq_min, freq_max, Nfreq), decimation_factor); } dft_force add_dft_force(const volume_list *where, const std::vector &freq, - int decimation_factor = 1) { + int decimation_factor = 0) { return add_dft_force(where, freq.data(), freq.size(), decimation_factor); } dft_force add_dft_force(const volume_list *where, const double *freq, size_t Nfreq, - int decimation_factor = 1); + int decimation_factor = 0); // near2far.cpp dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, - int Nfreq, int decimation_factor = 1, int Nperiods = 1) { + int Nfreq, int decimation_factor = 0, int Nperiods = 1) { return add_dft_near2far(where, linspace(freq_min, freq_max, Nfreq), decimation_factor, Nperiods); } dft_near2far add_dft_near2far(const volume_list *where, const std::vector &freq, - int decimation_factor = 1, int Nperiods = 1) { + int decimation_factor = 0, int Nperiods = 1) { return add_dft_near2far(where, freq.data(), freq.size(), decimation_factor, Nperiods); } dft_near2far add_dft_near2far(const volume_list *where, const double *freq, size_t Nfreq, - int decimation_factor = 1, int Nperiods = 1); + int decimation_factor = 0, int Nperiods = 1); // monitor.cpp std::complex get_chi1inv(component, direction, const vec &loc, double frequency = 0, bool parallel = true) const; diff --git a/src/sources.cpp b/src/sources.cpp index 65eda579c..3c814bd47 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -105,6 +105,13 @@ std::complex gaussian_src_time::fourier_transform(const double f) { return width * polar(1.0, omega * peak_time) * exp(-0.5 * delta * delta); } +// bandwidth (in frequency units, not angular frequency) of the +// continuous Fourier transform of the Gaussian source function +// when it has decayed by a tolerance tol below its peak value +double gaussian_src_time::get_fwidth(double tol) const { + return sqrt(-2.0 * log(tol)) / (width * pi); +} + bool gaussian_src_time::is_equal(const src_time &t) const { const gaussian_src_time *tp = dynamic_cast(&t); if (tp) diff --git a/tests/harmonics.cpp b/tests/harmonics.cpp index 81a8eb9ec..2f318c035 100644 --- a/tests/harmonics.cpp +++ b/tests/harmonics.cpp @@ -45,9 +45,12 @@ void harmonics(double freq, double chi2, double chi3, double J, double &A2, doub f.add_point_source(Ex, src, vec(-0.5 * sz + dpml), J); vec fpt(0.5 * sz - dpml - 0.5); - dft_flux d1 = f.add_dft_flux(Z, volume(fpt), freq, freq, 1); - dft_flux d2 = f.add_dft_flux(Z, volume(fpt), 2 * freq, 2 * freq, 1); - dft_flux d3 = f.add_dft_flux(Z, volume(fpt), 3 * freq, 3 * freq, 1); + dft_flux d1 = f.add_dft_flux(Z, volume(fpt), freq, freq, 1, true, true, + 1 /* decimation_factor */); + dft_flux d2 = f.add_dft_flux(Z, volume(fpt), 2 * freq, 2 * freq, 1, true, true, + 1 /* decimation_factor */); + dft_flux d3 = f.add_dft_flux(Z, volume(fpt), 3 * freq, 3 * freq, 1, true, true, + 1 /* decimation_factor */); double emax = 0;