From 3f3ccf5c506e133d74a17e4801d34347f89b2d62 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 13 Aug 2021 19:08:34 -0700 Subject: [PATCH 01/11] automatic DFT decimation --- doc/docs/Python_User_Interface.md | 75 ++++++++++++----------- python/adjoint/objective.py | 6 +- python/adjoint/optimization_problem.py | 2 +- python/simulation.py | 83 ++++++++++++++------------ src/dft.cpp | 15 +++++ src/meep.hpp | 52 ++++++++-------- tests/harmonics.cpp | 6 +- 7 files changed, 134 insertions(+), 105 deletions(-) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index c539ea90b..747952e60 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,12 +730,13 @@ 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 +`yee_grid` to `True`. To reduce the memory-bandwidth burden of accumulating +DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` +is 0 (the default), this value is automatically determined. It can be turned off by +setting it to 1. 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.
@@ -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,12 @@ 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. If `decimation_factor` +is 0 (the default), this value is automatically determined. It can be turned off by +setting it to 1. 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.
@@ -1361,7 +1364,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 +1494,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): ```
@@ -1502,12 +1505,13 @@ if they have not yet been initialized), telling Meep to accumulate the appropria 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 +below to get the energy spectrum, etcetera. To reduce the memory-bandwidth burden of accumulating +DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` +is 0 (the default), this value is automatically determined. It can be turned off by +setting it to 1. 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.
@@ -1723,7 +1727,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,11 +1739,12 @@ 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 +DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` +is 0 (the default), this value is automatically determined. It can be turned off by +setting it to 1. 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.
@@ -2008,7 +2013,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,11 +2025,12 @@ 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 +accumulating DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` +is 0 (the default), this value is automatically determined. It can be turned off by +setting it to 1. 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.
@@ -4266,8 +4272,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 framework 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/simulation.py b/python/simulation.py index 62a8b54a4..e5bd8b2f5 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,12 +2423,13 @@ 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 + `yee_grid` to `True`. To reduce the memory-bandwidth burden of accumulating + DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` + is 0 (the default), this value is automatically determined. It can be turned off by + setting it to 1. 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. """ components = args[0] @@ -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 + accumulating DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` + is 0 (the default), this value is automatically determined. It can be turned off by + setting it to 1. 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. """ 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,25 +2521,26 @@ 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 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 + below to get the energy spectrum, etcetera. To reduce the memory-bandwidth burden of accumulating + DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` + is 0 (the default), this value is automatically determined. It can be turned off by + setting it to 1. 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. """ 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 + DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` + is 0 (the default), this value is automatically determined. It can be turned off by + setting it to 1. 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. """ 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,18 @@ 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. If `decimation_factor` + is 0 (the default), this value is automatically determined. It can be turned off by + setting it to 1. 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. """ 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 +2881,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/src/dft.cpp b/src/dft.cpp index cf5177caa..e48012d6f 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -174,7 +174,22 @@ 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 src_freq_max = 0; + for (src_time *s = sources; s; s = s->next) { + if (s->frequency().real()+0.5*s->fwidth() > src_freq_max) + src_freq_max = s->frequency().real()+0.5*s->fwidth(); + } + if (src_freq_max > 0) + decimation_factor = 1/(dt*(freq[Nfreq-1] + src_freq_max)); + if (decimation_factor > 2) + decimation_factor -= 1; + 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..c3f002d7f 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 fwidth() 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 fwidth() const { return 1/width; }; 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 fwidth() const { return 1/width; }; 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 fwidth() 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/tests/harmonics.cpp b/tests/harmonics.cpp index 81a8eb9ec..945cd6b29 100644 --- a/tests/harmonics.cpp +++ b/tests/harmonics.cpp @@ -45,9 +45,9 @@ 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); + dft_flux d2 = f.add_dft_flux(Z, volume(fpt), 2 * freq, 2 * freq, 1, true, true, 1); + dft_flux d3 = f.add_dft_flux(Z, volume(fpt), 3 * freq, 3 * freq, 1, true, true, 1); double emax = 0; From 3e28428dc47b2e9dec417f7092501226fc32de0e Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sat, 14 Aug 2021 03:56:43 +0000 Subject: [PATCH 02/11] disable automatic decimation for certain tests --- python/tests/test_3rd_harm_1d.py | 6 +++--- python/tests/test_adjoint_jax.py | 3 ++- python/tests/test_bend_flux.py | 5 +++-- python/tests/test_cavity_farfield.py | 3 ++- python/tests/test_dft_energy.py | 6 ++++-- python/tests/test_dft_fields.py | 3 ++- python/tests/test_force.py | 2 +- python/tests/test_holey_wvg_cavity.py | 3 ++- python/tests/test_mode_decomposition.py | 3 ++- src/dft.cpp | 8 +++----- 10 files changed, 24 insertions(+), 18 deletions(-) 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_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_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_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/src/dft.cpp b/src/dft.cpp index e48012d6f..624de5ce8 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -181,12 +181,10 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, if (s->frequency().real()+0.5*s->fwidth() > src_freq_max) src_freq_max = s->frequency().real()+0.5*s->fwidth(); } - if (src_freq_max > 0) + if (src_freq_max > 0) { decimation_factor = 1/(dt*(freq[Nfreq-1] + src_freq_max)); - if (decimation_factor > 2) - decimation_factor -= 1; - else - decimation_factor = 1; + } + if (decimation_factor > 1) decimation_factor -= 1; } data.decimation_factor = decimation_factor; From 430487083a8aade302bea9f311f971ebbb9b0334 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sat, 14 Aug 2021 06:14:57 +0000 Subject: [PATCH 03/11] fixes --- doc/docs/Python_User_Interface.md | 2 +- python/adjoint/utils.py | 1 + python/geom.py | 2 +- python/tests/test_geom.py | 4 ++-- python/tests/test_visualization.py | 8 ++++---- src/dft.cpp | 9 +++++++-- src/meep.hpp | 8 ++++---- tests/harmonics.cpp | 9 ++++++--- 8 files changed, 26 insertions(+), 17 deletions(-) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 747952e60..7100b6bd8 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -4272,7 +4272,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. -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/adjoint/utils.py b/python/adjoint/utils.py index 110a4c4c2..be2abe4c5 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=1 ) 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/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_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 624de5ce8..2d2ccfacd 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -178,8 +178,13 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, if (decimation_factor == 0) { double src_freq_max = 0; for (src_time *s = sources; s; s = s->next) { - if (s->frequency().real()+0.5*s->fwidth() > src_freq_max) - src_freq_max = s->frequency().real()+0.5*s->fwidth(); + if (s->get_width() == 0) { + decimation_factor = 1; + } else { + if (s->frequency().real()+0.5/s->get_width() > src_freq_max) { + src_freq_max = s->frequency().real()+0.5/s->get_width(); + } + } } if (src_freq_max > 0) { decimation_factor = 1/(dt*(freq[Nfreq-1] + src_freq_max)); diff --git a/src/meep.hpp b/src/meep.hpp index c3f002d7f..36fe17d4b 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -973,7 +973,7 @@ class src_time { return 1; } virtual std::complex frequency() const { return 0.0; } - virtual double fwidth() const { return 0.0; } + virtual double get_width() const { return 0.0; } virtual void set_frequency(std::complex f) { (void)f; } private: @@ -995,7 +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 fwidth() const { return 1/width; }; + virtual double get_width() const { return width; }; virtual void set_frequency(std::complex f) { freq = real(f); } std::complex fourier_transform(const double f); @@ -1016,7 +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 fwidth() const { return 1/width; }; + virtual double get_width() const { return 0.0; }; virtual void set_frequency(std::complex f) { freq = f; } private: @@ -1049,7 +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 fwidth() const { return 0.0; }; + virtual double get_width() const { return 0.0; }; virtual void set_frequency(std::complex f) { freq = f; } private: diff --git a/tests/harmonics.cpp b/tests/harmonics.cpp index 945cd6b29..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, true, true, 1); - dft_flux d2 = f.add_dft_flux(Z, volume(fpt), 2 * freq, 2 * freq, 1, true, true, 1); - dft_flux d3 = f.add_dft_flux(Z, volume(fpt), 3 * freq, 3 * freq, 1, true, true, 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; From 89a31e86968353da089e2679df8eb54958b1999a Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sat, 14 Aug 2021 08:06:33 -0700 Subject: [PATCH 04/11] disable decimation for more tests --- python/tests/test_chunks.py | 4 ++-- python/tests/test_divide_mpi_processes.py | 5 ++--- python/tests/test_refl_angular.py | 4 ++-- 3 files changed, 6 insertions(+), 7 deletions(-) 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_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_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)) From 64592493b65535fec4304f8e9561da9ce3e98a87 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Tue, 17 Aug 2021 17:51:30 -0700 Subject: [PATCH 05/11] compute the maximum frequency of the DFT monitor rather than use the last array element --- doc/docs/Python_User_Interface.md | 52 ++++++++++++++++++------------- python/adjoint/utils.py | 2 +- python/simulation.py | 52 ++++++++++++++++++------------- src/dft.cpp | 10 ++++-- 4 files changed, 70 insertions(+), 46 deletions(-) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 7100b6bd8..a3bf7274f 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -731,10 +731,12 @@ default routine interpolates the Fourier-transformed fields at the center of eac 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` can be specified. If `decimation_factor` -is 0 (the default), this value is automatically determined. It can be turned off by -setting it to 1. DFT field values are updated every `decimation_factor` timesteps. -Use this feature with care, as the decimated timeseries may be corrupted by +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. @@ -1111,10 +1113,12 @@ 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` can be specified. If `decimation_factor` -is 0 (the default), this value is automatically determined. It can be turned off by -setting it to 1. DFT field values are updated every `decimation_factor` timesteps. -Use this feature with care, as the decimated timeseries may be corrupted by +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. @@ -1505,11 +1509,13 @@ if they have not yet been initialized), telling Meep to accumulate the appropria 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` can be specified. If `decimation_factor` -is 0 (the default), this value is automatically determined. It can be turned off by -setting it to 1. DFT field values are updated every `decimation_factor` timesteps. -Use this feature with care, as the decimated timeseries may be corrupted by +below to get the energy spectrum, etcetera. 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. 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. @@ -1739,10 +1745,12 @@ 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 -DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` -is 0 (the default), this value is automatically determined. It can be turned off by -setting it to 1. DFT field values are updated every `decimation_factor` timesteps. -Use this feature with care, as the decimated timeseries may be corrupted by +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. @@ -2025,10 +2033,12 @@ 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` can be specified. If `decimation_factor` -is 0 (the default), this value is automatically determined. It can be turned off by -setting it to 1. DFT field values are updated every `decimation_factor` timesteps. -Use this feature with care, as the decimated timeseries may be corrupted by +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. diff --git a/python/adjoint/utils.py b/python/adjoint/utils.py index be2abe4c5..decbe9e11 100644 --- a/python/adjoint/utils.py +++ b/python/adjoint/utils.py @@ -69,7 +69,7 @@ def install_design_region_monitors( frequencies, where=design_region.volume, yee_grid=True, - decimation_factor=1 + decimation_factor=0 ) for design_region in design_regions ] return design_region_monitors diff --git a/python/simulation.py b/python/simulation.py index e5bd8b2f5..89c788be2 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -2424,10 +2424,12 @@ def add_dft_fields(self, *args, **kwargs): 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` can be specified. If `decimation_factor` - is 0 (the default), this value is automatically determined. It can be turned off by - setting it to 1. DFT field values are updated every `decimation_factor` timesteps. - Use this feature with care, as the decimated timeseries may be corrupted by + 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. @@ -2496,10 +2498,12 @@ 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` can be specified. If `decimation_factor` - is 0 (the default), this value is automatically determined. It can be turned off by - setting it to 1. DFT field values are updated every `decimation_factor` timesteps. - Use this feature with care, as the decimated timeseries may be corrupted by + 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. @@ -2528,11 +2532,13 @@ def add_energy(self, *args, **kwargs): 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` can be specified. If `decimation_factor` - is 0 (the default), this value is automatically determined. It can be turned off by - setting it to 1. DFT field values are updated every `decimation_factor` timesteps. - Use this feature with care, as the decimated timeseries may be corrupted by + below to get the energy spectrum, etcetera. 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. 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. @@ -2751,10 +2757,12 @@ 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 - DFT fields, an integer `decimation_factor` can be specified. If `decimation_factor` - is 0 (the default), this value is automatically determined. It can be turned off by - setting it to 1. DFT field values are updated every `decimation_factor` timesteps. - Use this feature with care, as the decimated timeseries may be corrupted by + 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. @@ -2857,10 +2865,12 @@ 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` can be specified. If `decimation_factor` - is 0 (the default), this value is automatically determined. It can be turned off by - setting it to 1. DFT field values are updated every `decimation_factor` timesteps. - Use this feature with care, as the decimated timeseries may be corrupted by + 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. diff --git a/src/dft.cpp b/src/dft.cpp index 2d2ccfacd..00b420a7f 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -186,10 +186,14 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, } } } - if (src_freq_max > 0) { - decimation_factor = 1/(dt*(freq[Nfreq-1] + src_freq_max)); + double freq_max = 0; + for (size_t i = 0; i < Nfreq; ++i) { + if (freq[i] > freq_max) + freq_max = freq[i]; + } + if ((freq_max > 0) && (src_freq_max > 0)) { + decimation_factor = 1/(dt*(freq_max + src_freq_max)); } - if (decimation_factor > 1) decimation_factor -= 1; } data.decimation_factor = decimation_factor; From 90f21e14f8cd97687379b38d0d209763fb91dad0 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Wed, 18 Aug 2021 12:51:28 -0700 Subject: [PATCH 06/11] replace get_width with get_fwidth and use analytic formula for gaussian_src_time --- doc/docs/Python_User_Interface.md | 16 ++++------------ python/simulation.py | 16 ++++------------ src/dft.cpp | 7 ++++--- src/meep.hpp | 8 ++++---- src/sources.cpp | 6 ++++++ 5 files changed, 22 insertions(+), 31 deletions(-) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index a3bf7274f..33f59f655 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -737,9 +737,7 @@ 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. +[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies.
@@ -1516,9 +1514,7 @@ 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. +[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. @@ -1751,9 +1747,7 @@ 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. +[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. @@ -2039,9 +2033,7 @@ 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. +[aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. diff --git a/python/simulation.py b/python/simulation.py index 89c788be2..34231ab6e 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -2430,9 +2430,7 @@ def add_dft_fields(self, *args, **kwargs): [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. + [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. """ components = args[0] args = fix_dft_args(args, 1) @@ -2504,9 +2502,7 @@ def add_near2far(self, *args, **kwargs): [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. + [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. """ args = fix_dft_args(args, 0) freq = args[0] @@ -2539,9 +2535,7 @@ def add_energy(self, *args, **kwargs): [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. + [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. """ args = fix_dft_args(args, 0) freq = args[0] @@ -2763,9 +2757,7 @@ def add_force(self, *args, **kwargs): [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. + [aliasing](https://en.wikipedia.org/wiki/Aliasing) of high frequencies. """ args = fix_dft_args(args, 0) freq = args[0] diff --git a/src/dft.cpp b/src/dft.cpp index 00b420a7f..80498ce24 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -176,13 +176,14 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, 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_width() == 0) { + if (s->get_fwidth(tol) == 0) { decimation_factor = 1; } else { - if (s->frequency().real()+0.5/s->get_width() > src_freq_max) { - src_freq_max = s->frequency().real()+0.5/s->get_width(); + if (s->frequency().real()+0.5*s->get_fwidth(tol) > src_freq_max) { + src_freq_max = s->frequency().real()+0.5*s->get_fwidth(tol); } } } diff --git a/src/meep.hpp b/src/meep.hpp index 36fe17d4b..7dc9fdf93 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -973,7 +973,7 @@ class src_time { return 1; } virtual std::complex frequency() const { return 0.0; } - virtual double get_width() const { return 0.0; } + virtual double get_fwidth(double tol) const { return 0.0; } virtual void set_frequency(std::complex f) { (void)f; } private: @@ -995,7 +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_width() const { return width; }; + virtual double get_fwidth(double tol) const; virtual void set_frequency(std::complex f) { freq = real(f); } std::complex fourier_transform(const double f); @@ -1016,7 +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_width() const { return 0.0; }; + virtual double get_fwidth(double tol) const { return 0.0; }; virtual void set_frequency(std::complex f) { freq = f; } private: @@ -1049,7 +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_width() const { return 0.0; }; + virtual double get_fwidth(double tol) const { return 0.0; }; virtual void set_frequency(std::complex f) { freq = f; } private: diff --git a/src/sources.cpp b/src/sources.cpp index 65eda579c..7397ec052 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -105,6 +105,12 @@ std::complex gaussian_src_time::fourier_transform(const double f) { return width * polar(1.0, omega * peak_time) * exp(-0.5 * delta * delta); } +// bandwidth of the Fourier transform of the Gaussian source function +// when it has decayed by tolerance tol below its peak value +double gaussian_src_time::get_fwidth(double tol) const { + return 2/width * sqrt(-2.0 * log(tol)); +} + bool gaussian_src_time::is_equal(const src_time &t) const { const gaussian_src_time *tp = dynamic_cast(&t); if (tp) From 8fb8b061c2c48476960e3dee1f1684e93e93fccd Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 19 Aug 2021 10:18:33 -0700 Subject: [PATCH 07/11] add missing 2pi factor to gaussian_src_time::get_fwidth --- src/sources.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/sources.cpp b/src/sources.cpp index 7397ec052..a0049834e 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -105,10 +105,11 @@ std::complex gaussian_src_time::fourier_transform(const double f) { return width * polar(1.0, omega * peak_time) * exp(-0.5 * delta * delta); } -// bandwidth of the Fourier transform of the Gaussian source function -// when it has decayed by tolerance tol below its peak value +// bandwidth (in frequency units, not angular frequency) of the +// 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 2/width * sqrt(-2.0 * log(tol)); + return sqrt(-2.0 * log(tol)) / (width * 2 * pi); } bool gaussian_src_time::is_equal(const src_time &t) const { From 400fa18de3393e7b965f651443b3db92bf23c76b Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 19 Aug 2021 12:42:41 -0700 Subject: [PATCH 08/11] add factor of 2 to gaussian_src_time::get_fwidth(double tol) --- src/sources.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/sources.cpp b/src/sources.cpp index a0049834e..3c814bd47 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -106,10 +106,10 @@ std::complex gaussian_src_time::fourier_transform(const double f) { } // bandwidth (in frequency units, not angular frequency) of the -// Fourier transform of the Gaussian source function when it has -// decayed by a tolerance tol below its peak value +// 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 * 2 * pi); + return sqrt(-2.0 * log(tol)) / (width * pi); } bool gaussian_src_time::is_equal(const src_time &t) const { From c8af01d85619fccb52462dc6d59d3a41b5ccb537 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 19 Aug 2021 19:03:06 -0700 Subject: [PATCH 09/11] various fixes --- src/dft.cpp | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index 80498ce24..dbc356544 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -179,22 +179,18 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, 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 { - if (s->frequency().real()+0.5*s->get_fwidth(tol) > src_freq_max) { - src_freq_max = s->frequency().real()+0.5*s->get_fwidth(tol); - } - } + 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) { - if (freq[i] > freq_max) - freq_max = freq[i]; - } - if ((freq_max > 0) && (src_freq_max > 0)) { - decimation_factor = 1/(dt*(freq_max + src_freq_max)); - } + 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 = int(std::ceil(1/(dt*(freq_max + src_freq_max)))); + else + decimation_factor = 1; } data.decimation_factor = decimation_factor; From e4815421f0e020c7982b3cab7f3069bd436db2dd Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 20 Aug 2021 02:49:43 +0000 Subject: [PATCH 10/11] disable automatic decimation for test_n2f_periodic.py --- python/tests/test_n2f_periodic.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/python/tests/test_n2f_periodic.py b/python/tests/test_n2f_periodic.py index 5cac78d5c..840cdb990 100644 --- a/python/tests/test_n2f_periodic.py +++ b/python/tests/test_n2f_periodic.py @@ -49,8 +49,10 @@ def test_nea2far_periodic(self): sources=sources, symmetries=symmetries) - n2f_obj = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy)), nperiods=10) - dft_obj = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=dft_pt, size=mp.Vector3(y=sy)) + n2f_obj = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy)), + nperiods=10, decimation_factor=1) + dft_obj = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=dft_pt, size=mp.Vector3(y=sy), + decimation_factor=1) sim.run(until_after_sources=300) From 898d438505050fc343384e1eee10f4994ecac817 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Mon, 23 Aug 2021 15:41:35 -0700 Subject: [PATCH 11/11] switch ceil to floor when setting decimation_factor --- python/tests/test_n2f_periodic.py | 6 ++---- src/dft.cpp | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/python/tests/test_n2f_periodic.py b/python/tests/test_n2f_periodic.py index 840cdb990..5cac78d5c 100644 --- a/python/tests/test_n2f_periodic.py +++ b/python/tests/test_n2f_periodic.py @@ -49,10 +49,8 @@ def test_nea2far_periodic(self): sources=sources, symmetries=symmetries) - n2f_obj = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy)), - nperiods=10, decimation_factor=1) - dft_obj = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=dft_pt, size=mp.Vector3(y=sy), - decimation_factor=1) + n2f_obj = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy)), nperiods=10) + dft_obj = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=dft_pt, size=mp.Vector3(y=sy)) sim.run(until_after_sources=300) diff --git a/src/dft.cpp b/src/dft.cpp index dbc356544..0ecf4825f 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -188,7 +188,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, 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 = int(std::ceil(1/(dt*(freq_max + src_freq_max)))); + decimation_factor = std::max(1, int(std::floor(1/(dt*(freq_max + src_freq_max))))); else decimation_factor = 1; }