From 8f8599347b91537490d9ea0bf60f95fafdf39884 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 6 Aug 2021 07:09:03 -0700 Subject: [PATCH] decimation option for remaining DFT functions (#1720) * add decimation factor to add_dft_flux * some fixes * fixes * fixes, unit tests, docs * update markdown pages for user manual * increase tolerance of failing unit test in test_bend_flux.py * fix add_dft_near2far function header --- doc/docs/Python_User_Interface.md | 107 +++++++++++++++++++----- python/simulation.py | 106 +++++++++++++++-------- python/tests/test_3rd_harm_1d.py | 2 - python/tests/test_antenna_radiation.py | 14 ++-- python/tests/test_bend_flux.py | 40 ++++++--- python/tests/test_dft_energy.py | 10 ++- python/tests/test_dft_fields.py | 15 ++-- python/tests/test_force.py | 3 +- python/tests/test_mode_decomposition.py | 13 +++ python/tests/test_n2f_periodic.py | 6 +- src/dft.cpp | 35 +++++--- src/meep.hpp | 82 +++++++++++------- src/near2far.cpp | 5 +- src/stress.cpp | 17 ++-- 14 files changed, 311 insertions(+), 144 deletions(-) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index ad351dc69..72346478c 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -534,10 +534,29 @@ def mean_time_spent_on(self, time_sink):
Return the mean time spent by all processes for a type of work `time_sink` which -can be one of ten integer values `0`-`9`: (`0`) connecting chunks, (`1`) time stepping, -(`2`) copying boundaries, (`3`) MPI all-to-all communication/synchronization, -(`4`) MPI one-to-one communication, (`5`) field output, (`6`) Fourier transforming, -(`7`) MPB mode solver, (`8`) near-to-far field transformation, and (`9`) other. +can be one of the following integer constants: +* meep.Stepping ("time stepping") +* meep.Connecting ("connecting chunks") +* meep.Boundaries ("copying boundaries") +* meep.MpiAllTime ("all-all communication") +* meep.MpiOneTime ("1-1 communication") +* meep.FieldOutput ("outputting fields") +* meep.FourierTransforming ("Fourier transforming") +* meep.MPBTime ("MPB mode solver") +* meep.GetFarfieldsTime ("far-field transform") +* meep.FieldUpdateB ("updating B field") +* meep.FieldUpdateH ("updating H field") +* meep.FieldUpdateD ("updating D field") +* meep.FieldUpdateE ("updating E field") +* meep.BoundarySteppingB ("boundary stepping B") +* meep.BoundarySteppingWH ("boundary stepping WH") +* meep.BoundarySteppingPH ("boundary stepping PH") +* meep.BoundarySteppingH ("boundary stepping H") +* meep.BoundarySteppingD ("boundary stepping D") +* meep.BoundarySteppingWE ("boundary stepping WE") +* meep.BoundarySteppingPE ("boundary stepping PE") +* meep.BoundarySteppingE ("boundary stepping E") +* meep.Other ("everything else")
@@ -726,16 +745,16 @@ fields 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 over the `Volume` specified by `where` (default to the entire cell). 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 +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 +`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 experimental feature with care, as the decimated -timeseries may be corruped by aliasing of high frequencies. The choice -of decimation factor should take the properties of all sources in the -simulation and the frequency range of the DFT field monitor into account. +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. @@ -1009,7 +1028,9 @@ def reset_meep(self):
Reset all of Meep's parameters, deleting the fields, structures, etcetera, from -memory as if you had not run any computations. +memory as if you had not run any computations. If the `num_chunks` or `chunk_layout` +attributes have been modified internally, they are reset to their original +values passed in at instantiation.
@@ -1094,8 +1115,8 @@ Given a bunch of [`FluxRegion`](#fluxregion) objects, you can tell Meep to accum
```python -def add_flux(self, *args): -def add_flux(fcen, df, nfreq, freq, FluxRegions...): +def add_flux(self, *args, **kwargs): +def add_flux(fcen, df, nfreq, freq, FluxRegions, decimation_factor=1): ```
@@ -1106,7 +1127,13 @@ 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 a *flux object*, which you can pass to the functions -below to get the flux spectrum, etcetera. +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 +in the simulation as well as the frequency range of the DFT field monitor.
@@ -1482,8 +1509,8 @@ The usage is similar to the flux spectra: you define a set of [`EnergyRegion`](#
```python -def add_energy(self, *args): -def add_energy(fcen, df, nfreq, freq, EnergyRegions...): +def add_energy(self, *args, **kwargs): +def add_energy(fcen, df, nfreq, freq, EnergyRegions, decimation_factor=1): ```
@@ -1494,7 +1521,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. +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.
@@ -1708,8 +1741,8 @@ The usage is similar to the [flux spectra](Python_Tutorials/Basics.md#transmitta
```python -def add_force(self, *args): -def add_force(fcen, df, nfreq, freq, ForceRegions...): +def add_force(self, *args, **kwargs): +def add_force(fcen, df, nfreq, freq, ForceRegions, decimation_factor=1): ```
@@ -1720,7 +1753,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 a `force`object, which you can pass to the functions -below to get the force spectrum, etcetera. +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.
@@ -1988,7 +2027,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): +def add_near2far(fcen, df, nfreq, freq, Near2FarRegions, nperiods=1, decimation_factor=1): ```
@@ -1999,7 +2038,13 @@ 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 a `near2far` object, which you can pass -to the functions below to get the far fields. +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.
@@ -7366,6 +7411,22 @@ of processes.
+ + +
+ +```python +def print(self): +``` + +
+ +Pretty-prints the tree structure of the BinaryPartition object. + +
+ +
+ Miscellaneous Functions Reference --------------------------------- diff --git a/python/simulation.py b/python/simulation.py index f7fe7effd..9eca33bda 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -2418,16 +2418,16 @@ def add_dft_fields(self, *args, **kwargs): `fcen-df/2` to `fcen+df/2` or an array/list `freq` for arbitrarily spaced frequencies over the `Volume` specified by `where` (default to the entire cell). 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 + 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 + `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 experimental feature with care, as the decimated - timeseries may be corruped by aliasing of high frequencies. The choice - of decimation factor should take the properties of all sources in the - simulation and the frequency range of the DFT field monitor into account. + 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] args = fix_dft_args(args, 1) @@ -2485,50 +2485,66 @@ def get_dft_data(self, dft_chunk): def add_near2far(self, *args, **kwargs): """ - `add_near2far(fcen, df, nfreq, freq, Near2FarRegions..., nperiods=1)` ##sig + `add_near2far(fcen, df, nfreq, freq, Near2FarRegions, nperiods=1, decimation_factor=1)` ##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 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 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. """ args = fix_dft_args(args, 0) freq = args[0] near2fars = args[1:] nperiods = kwargs.get('nperiods', 1) - n2f = DftNear2Far(self._add_near2far, [freq, nperiods, near2fars]) + decimation_factor = kwargs.get('decimation_factor', 1) + n2f = DftNear2Far(self._add_near2far, [freq, nperiods, near2fars, decimation_factor]) self.dft_objects.append(n2f) return n2f - def _add_near2far(self, freq, nperiods, near2fars): + def _add_near2far(self, freq, nperiods, near2fars, decimation_factor): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_near2far, freq, near2fars, nperiods) + return self._add_fluxish_stuff(self.fields.add_dft_near2far, freq, near2fars, + decimation_factor, nperiods) - def add_energy(self, *args): + def add_energy(self, *args, **kwargs): """ - `add_energy(fcen, df, nfreq, freq, EnergyRegions...)` ##sig + `add_energy(fcen, df, nfreq, freq, EnergyRegions, decimation_factor=1)` ##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. + 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. """ args = fix_dft_args(args, 0) freq = args[0] energys = args[1:] - en = DftEnergy(self._add_energy, [freq, energys]) + decimation_factor = kwargs.get('decimation_factor', 1) + en = DftEnergy(self._add_energy, [freq, energys, decimation_factor]) self.dft_objects.append(en) return en - def _add_energy(self, freq, energys): + def _add_energy(self, freq, energys, decimation_factor): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_energy, freq, energys) + return self._add_fluxish_stuff(self.fields.add_dft_energy, freq, energys, + decimation_factor) def _display_energy(self, name, func, energys): if energys: @@ -2720,28 +2736,36 @@ def load_minus_near2far_data(self, near2far, n2fdata): self.load_near2far_data(near2far, n2fdata) near2far.scale_dfts(complex(-1.0)) - def add_force(self, *args): + def add_force(self, *args, **kwargs): """ - `add_force(fcen, df, nfreq, freq, ForceRegions...)` ##sig + `add_force(fcen, df, nfreq, freq, ForceRegions, decimation_factor=1)` ##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 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. + 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. """ args = fix_dft_args(args, 0) freq = args[0] forces = args[1:] - force = DftForce(self._add_force, [freq, forces]) + decimation_factor = kwargs.get('decimation_factor', 1) + force = DftForce(self._add_force, [freq, forces, decimation_factor]) self.dft_objects.append(force) return force - def _add_force(self, freq, forces): + def _add_force(self, freq, forces, decimation_factor): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_force, freq, forces) + return self._add_fluxish_stuff(self.fields.add_dft_force, freq, forces, + decimation_factor) def display_forces(self, *forces): """ @@ -2817,28 +2841,36 @@ def load_minus_force_data(self, force, fdata): self.load_force_data(force, fdata) force.scale_dfts(complex(-1.0)) - def add_flux(self, *args): + def add_flux(self, *args, **kwargs): """ - `add_flux(fcen, df, nfreq, freq, FluxRegions...)` ##sig + `add_flux(fcen, df, nfreq, freq, FluxRegions, decimation_factor=1)` ##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 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. + 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 + 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:] - flux = DftFlux(self._add_flux, [freq, fluxes]) + decimation_factor = kwargs.get('decimation_factor', 1) + flux = DftFlux(self._add_flux, [freq, fluxes, decimation_factor]) self.dft_objects.append(flux) return flux - def _add_flux(self, freq, fluxes): + def _add_flux(self, freq, fluxes, decimation_factor): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_flux, freq, fluxes) + return self._add_fluxish_stuff(self.fields.add_dft_flux, + freq, fluxes, decimation_factor) def add_mode_monitor(self, *args, **kwargs): """ @@ -2849,12 +2881,13 @@ def add_mode_monitor(self, *args, **kwargs): args = fix_dft_args(args, 0) freq = args[0] fluxes = args[1:] + decimation_factor = kwargs.get('decimation_factor', 1) yee_grid = kwargs.get("yee_grid", False) - flux = DftFlux(self._add_mode_monitor, [freq, fluxes, yee_grid]) + flux = DftFlux(self._add_mode_monitor, [freq, fluxes, yee_grid, decimation_factor]) self.dft_objects.append(flux) return flux - def _add_mode_monitor(self, freq, fluxes, yee_grid): + def _add_mode_monitor(self, freq, fluxes, yee_grid, decimation_factor): if self.fields is None: self.init_sim() @@ -2867,7 +2900,7 @@ def _add_mode_monitor(self, freq, fluxes, yee_grid): d0 = region.direction d = self.fields.normal_direction(v.swigobj) if d0 < 0 else d0 - return self.fields.add_mode_monitor(d, v.swigobj, freq, centered_grid) + return self.fields.add_mode_monitor(d, v.swigobj, freq, centered_grid, decimation_factor) def display_fluxes(self, *fluxes): """ @@ -3082,7 +3115,7 @@ def solve_eigfreq(self, tol=1e-7, maxiters=100, self.fields.solve_cw(cwtol, cwmaxiters, guessfreq, L, eigfreq, tol, maxiters) return eigfreq.item() - def _add_fluxish_stuff(self, add_dft_stuff, freq, stufflist, *args): + def _add_fluxish_stuff(self, add_dft_stuff, freq, stufflist, decimation_factor, *args): vol_list = None for s in stufflist: @@ -3094,8 +3127,7 @@ def _add_fluxish_stuff(self, add_dft_stuff, freq, stufflist, *args): v2 = Volume(center=s.center, size=s.size, dims=self.dimensions, is_cylindrical=self.is_cylindrical).swigobj vol_list = mp.make_volume_list(v2, c, s.weight, vol_list) - - stuff = add_dft_stuff(vol_list, freq, *args) + stuff = add_dft_stuff(vol_list, freq, decimation_factor, *args) vol_list.__swig_destroy__(vol_list) return stuff @@ -3618,7 +3650,7 @@ def change_sources(self, new_sources): def reset_meep(self): """ Reset all of Meep's parameters, deleting the fields, structures, etcetera, from - memory as if you had not run any computations. If the num_chunks or chunk_layout + memory as if you had not run any computations. If the `num_chunks` or `chunk_layout` attributes have been modified internally, they are reset to their original values passed in at instantiation. """ diff --git a/python/tests/test_3rd_harm_1d.py b/python/tests/test_3rd_harm_1d.py index 6c0b52b05..f0bf9269f 100644 --- a/python/tests/test_3rd_harm_1d.py +++ b/python/tests/test_3rd_harm_1d.py @@ -1,5 +1,3 @@ -from __future__ import division - import unittest import meep as mp from utils import ApproxComparisonTestCase diff --git a/python/tests/test_antenna_radiation.py b/python/tests/test_antenna_radiation.py index df64f90f6..9aec952de 100644 --- a/python/tests/test_antenna_radiation.py +++ b/python/tests/test_antenna_radiation.py @@ -41,7 +41,8 @@ def test_farfield(self): mp.Near2FarRegion(mp.Vector3(y=0.5*sxy), size=mp.Vector3(sxy)), mp.Near2FarRegion(mp.Vector3(y=-0.5*sxy), size=mp.Vector3(sxy), weight=-1), mp.Near2FarRegion(mp.Vector3(0.5*sxy), size=mp.Vector3(y=sxy)), - mp.Near2FarRegion(mp.Vector3(-0.5*sxy), size=mp.Vector3(y=sxy), weight=-1)) + mp.Near2FarRegion(mp.Vector3(-0.5*sxy), size=mp.Vector3(y=sxy), weight=-1), + decimation_factor=20) flux_box = sim.add_flux(fcen, 0, 1, mp.FluxRegion(mp.Vector3(y=0.5*sxy), size=mp.Vector3(sxy)), @@ -69,11 +70,12 @@ def test_farfield(self): Pr = np.sqrt(np.square(Px)+np.square(Py)) far_flux_circle = np.sum(Pr)*2*np.pi*r/len(Pr) - rr = 20/fcen # length of far field square box - far_flux_square = (nearfield_box.flux(mp.Y, mp.Volume(center=mp.Vector3(y=0.5*rr), size=mp.Vector3(rr)), resolution)[0] - - nearfield_box.flux(mp.Y, mp.Volume(center=mp.Vector3(y=-0.5*rr), size=mp.Vector3(rr)), resolution)[0] - + nearfield_box.flux(mp.X, mp.Volume(center=mp.Vector3(0.5*rr), size=mp.Vector3(y=rr)), resolution)[0] - - nearfield_box.flux(mp.X, mp.Volume(center=mp.Vector3(-0.5*rr), size=mp.Vector3(y=rr)), resolution)[0]) + rr = 20/fcen # length of far-field square box + res_far = 20 # resolution of far-field square box + far_flux_square = (nearfield_box.flux(mp.Y, mp.Volume(center=mp.Vector3(y=0.5*rr), size=mp.Vector3(rr)), res_far)[0] + - nearfield_box.flux(mp.Y, mp.Volume(center=mp.Vector3(y=-0.5*rr), size=mp.Vector3(rr)), res_far)[0] + + nearfield_box.flux(mp.X, mp.Volume(center=mp.Vector3(0.5*rr), size=mp.Vector3(y=rr)), res_far)[0] + - nearfield_box.flux(mp.X, mp.Volume(center=mp.Vector3(-0.5*rr), size=mp.Vector3(y=rr)), res_far)[0]) print("flux:, {:.6f}, {:.6f}, {:.6f}".format(near_flux,far_flux_circle,far_flux_square)) diff --git a/python/tests/test_bend_flux.py b/python/tests/test_bend_flux.py index 70c05769a..17dff7950 100644 --- a/python/tests/test_bend_flux.py +++ b/python/tests/test_bend_flux.py @@ -1,13 +1,11 @@ -from __future__ import division - import os import unittest import numpy as np import meep as mp -from utils import compare_arrays +from utils import ApproxComparisonTestCase -class TestBendFlux(unittest.TestCase): +class TestBendFlux(ApproxComparisonTestCase): def init(self, no_bend=False, gdsii=False): sx = 16 @@ -65,10 +63,13 @@ def init(self, no_bend=False, gdsii=False): 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_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_decimated = self.sim.add_flux(np.linspace(fcen-0.5*df,fcen+0.5*df,nfreq), refl_fr, + decimation_factor=10) if no_bend: self.pt = mp.Vector3((sx / 2) - 1.5, wvg_ycen) @@ -81,6 +82,7 @@ def run_bend_flux(self, from_gdsii_file): self.sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, self.pt, 1e-3)) # Save flux data for use in real run below fdata = self.sim.get_flux_data(self.refl) + fdata_decimated = self.sim.get_flux_data(self.refl_decimated) expected = [ (0.1, 3.65231563251e-05, 3.68932495077e-05), @@ -105,16 +107,24 @@ def run_bend_flux(self, from_gdsii_file): (0.119191919192, 0.0254987474079, 0.0252348211592), ] - res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) + res = list(zip(mp.get_flux_freqs(self.trans), + mp.get_fluxes(self.trans), + mp.get_fluxes(self.refl))) - tolerance = 1e-2 if from_gdsii_file else 1e-3 - compare_arrays(self, np.array(expected), np.array(res[:20]), tol=tolerance) + res_decimated = list(zip(mp.get_flux_freqs(self.trans_decimated), + mp.get_fluxes(self.trans_decimated), + mp.get_fluxes(self.refl_decimated))) + + tol = 1e-6 if mp.is_single_precision() else 1e-8 + self.assertClose(np.array(expected), np.array(res[:20]), epsilon=tol) + self.assertClose(np.array(expected), np.array(res_decimated[:20]), epsilon=tol) # Real run self.sim = None self.init(gdsii=from_gdsii_file) # Load flux data obtained from normalization run self.sim.load_minus_flux_data(self.refl, fdata) + self.sim.load_minus_flux_data(self.refl_decimated, fdata_decimated) self.sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, self.pt, 1e-3)) expected = [ @@ -138,13 +148,19 @@ def run_bend_flux(self, from_gdsii_file): (0.11717171717171727, 0.005243320443599316, -0.003215529845495731), (0.11818181818181829, 0.0067654019326068, -0.004266367104375331), (0.11919191919191931, 0.008646855439680507, -0.005614491919262783), - ] - res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) + res = list(zip(mp.get_flux_freqs(self.trans), + mp.get_fluxes(self.trans), + mp.get_fluxes(self.refl))) + + res_decimated = list(zip(mp.get_flux_freqs(self.trans_decimated), + mp.get_fluxes(self.trans_decimated), + mp.get_fluxes(self.refl_decimated))) - tolerance = 1e-2 if from_gdsii_file else 1e-3 - compare_arrays(self, np.array(expected), np.array(res[:20]), tol=tolerance) + tol = 1e-3 + self.assertClose(np.array(expected), np.array(res[:20]), epsilon=tol) + self.assertClose(np.array(expected), np.array(res_decimated[:20]), epsilon=tol) def test_bend_flux(self): self.run_bend_flux(False) diff --git a/python/tests/test_dft_energy.py b/python/tests/test_dft_energy.py index 94bf01e53..856697a81 100644 --- a/python/tests/test_dft_energy.py +++ b/python/tests/test_dft_energy.py @@ -1,5 +1,3 @@ -from __future__ import division - import unittest import meep as mp @@ -26,6 +24,9 @@ def test_dft_energy(self): 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))) + energy_decimated = sim.add_energy(fsrc, 0, 1, + mp.EnergyRegion(center=mp.Vector3(3), size=mp.Vector3(y=5)), + decimation_factor=10) sim.run(until_after_sources=100) res = sim.get_eigenmode_coefficients(flux, [1], eig_parity=mp.ODD_Z+mp.EVEN_Y) @@ -39,6 +40,11 @@ def test_dft_energy(self): self.assertAlmostEqual(m_energy + e_energy, t_energy) self.assertAlmostEqual(ratio_vg, mode_vg, places=3) + e_energy_decimated = mp.get_electric_energy(energy_decimated)[0] + m_energy_decimated = mp.get_magnetic_energy(energy_decimated)[0] + self.assertAlmostEqual(e_energy, e_energy_decimated, places=1) + self.assertAlmostEqual(m_energy, m_energy_decimated, places=1) + if __name__ == '__main__': unittest.main() diff --git a/python/tests/test_dft_fields.py b/python/tests/test_dft_fields.py index fd11e3799..a66954abb 100644 --- a/python/tests/test_dft_fields.py +++ b/python/tests/test_dft_fields.py @@ -2,9 +2,10 @@ import h5py import numpy as np import meep as mp +from utils import ApproxComparisonTestCase import os -class TestDFTFields(unittest.TestCase): +class TestDFTFields(ApproxComparisonTestCase): @classmethod def setUpClass(cls): @@ -79,8 +80,9 @@ def test_get_dft_array(self): with h5py.File(os.path.join(self.temp_dir, 'thin-y-flux.h5'), 'r') as thin_y: thin_y_h5 = mp.complexarray(thin_y['ez_0.r'][()], thin_y['ez_0.i'][()]) - np.testing.assert_allclose(thin_x_array, thin_x_h5) - np.testing.assert_allclose(thin_y_array, thin_y_h5) + tol = 1e-6 + self.assertClose(thin_x_array, thin_x_h5, epsilon=tol) + self.assertClose(thin_y_array, thin_y_h5, epsilon=tol) # compare array data to HDF5 file content for fields and flux fields_arr = sim.get_dft_array(dft_fields, mp.Ez, 0) @@ -93,8 +95,9 @@ def test_get_dft_array(self): exp_fields = mp.complexarray(fields['ez_0.r'][()], fields['ez_0.i'][()]) exp_flux = mp.complexarray(flux['ez_0.r'][()], flux['ez_0.i'][()]) - np.testing.assert_allclose(exp_fields, fields_arr) - np.testing.assert_allclose(exp_flux, flux_arr) + tol = 1e-6 + self.assertClose(exp_fields, fields_arr, epsilon=tol) + self.assertClose(exp_flux, flux_arr, epsilon=tol) def test_decimated_dft_fields_are_almost_equal_to_undecimated_fields(self): sim = self.init() @@ -110,7 +113,7 @@ def test_decimated_dft_fields_are_almost_equal_to_undecimated_fields(self): expected_dft = sim.get_dft_array(undecimated_field, mp.Ez, 0) actual_dft = sim.get_dft_array(decimated_field, mp.Ez, 0) - np.testing.assert_allclose(expected_dft, actual_dft, rtol=5.e-2) + self.assertClose(expected_dft, actual_dft, epsilon=1e-3) if __name__ == '__main__': diff --git a/python/tests/test_force.py b/python/tests/test_force.py index b0a1dabb6..136202d0d 100644 --- a/python/tests/test_force.py +++ b/python/tests/test_force.py @@ -1,5 +1,4 @@ import unittest - import meep as mp @@ -22,6 +21,7 @@ def setUp(self): 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_decimated = self.sim.add_force(fcen, 0, 1, fr, decimation_factor=10) def test_force(self): @@ -35,6 +35,7 @@ def test_force(self): f = mp.get_forces(self.myforce) self.assertAlmostEqual(f[0], -0.11039089113393187) + self.assertAlmostEqual(f[0], mp.get_forces(self.myforce_decimated)[0]) if __name__ == '__main__': diff --git a/python/tests/test_mode_decomposition.py b/python/tests/test_mode_decomposition.py index 5d861fddb..b6abcd0e4 100644 --- a/python/tests/test_mode_decomposition.py +++ b/python/tests/test_mode_decomposition.py @@ -117,6 +117,10 @@ 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))) + 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)), + decimation_factor=10) sim.run(until_after_sources=30) @@ -124,10 +128,19 @@ def test_oblique_waveguide_backward_mode(self): coeff = sim.get_eigenmode_coefficients(mode,[1], direction=mp.NO_DIRECTION, kpoint_func=lambda f,n: kpoint).alpha[0,0,0] + flux_decimated = mp.get_fluxes(mode_decimated)[0] + coeff_decimated = sim.get_eigenmode_coefficients(mode_decimated,[1], + direction=mp.NO_DIRECTION, + kpoint_func=lambda f,n: kpoint).alpha[0,0,0] print("oblique-waveguide-flux:, {:.6f}, {:.6f}".format(-flux, abs(coeff)**2)) + print("oblique-waveguide-flux (decimated):, {:.6f}, {:.6f}".format(-flux_decimated, + abs(coeff_decimated)**2)) ## the magnitude of |flux| is 100.008731 and so we check two significant digits of accuracy self.assertAlmostEqual(-1,abs(coeff)**2/flux,places=2) + self.assertAlmostEqual(flux,flux_decimated,places=3) + self.assertAlmostEqual(coeff,coeff_decimated,places=3) + if __name__ == '__main__': unittest.main() diff --git a/python/tests/test_n2f_periodic.py b/python/tests/test_n2f_periodic.py index c9171369f..5cac78d5c 100644 --- a/python/tests/test_n2f_periodic.py +++ b/python/tests/test_n2f_periodic.py @@ -40,7 +40,7 @@ def test_nea2far_periodic(self): res = [20,25,30] norm = np.empty(3) - for j in range(3): + for j in range(3): sim = mp.Simulation(resolution=res[j], cell_size=mp.Vector3(sx,sy), boundary_layers=pml_layers, @@ -59,10 +59,10 @@ def test_nea2far_periodic(self): norm[j] = LA.norm(n2f_Ez['Ez']-dft_Ez[1:-1]) print("norm:, {}, {:.5f}".format(res[j],norm[j])) - sim.reset_meep() + sim.reset_meep() self.assertGreater(norm[0],norm[1]) self.assertGreater(norm[1],norm[2]) if __name__ == '__main__': - unittest.main() + unittest.main() diff --git a/src/dft.cpp b/src/dft.cpp index 8be4173eb..cf5177caa 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -158,7 +158,8 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, ive dft_chunk *fields::add_dft(component c, const volume &where, const double *freq, size_t Nfreq, bool include_dV_and_interp_weights, complex stored_weight, dft_chunk *chunk_next, bool sqrt_dV_and_interp_weights, - complex extra_weight, bool use_centered_grid, int vc, int decimation_factor) { + complex extra_weight, bool use_centered_grid, + int vc, int decimation_factor) { if (coordinate_mismatch(gv.dim, c)) return NULL; /* If you call add_dft before adding sources, it will do nothing @@ -423,7 +424,7 @@ void dft_flux::scale_dfts(complex scale) { } dft_flux fields::add_dft_flux(const volume_list *where_, const double *freq, size_t Nfreq, - bool use_symmetry, bool centered_grid) { + bool use_symmetry, bool centered_grid, int decimation_factor) { if (!where_) // handle empty list of volumes return dft_flux(Ex, Hy, NULL, NULL, freq, Nfreq, v, NO_DIRECTION, use_symmetry); @@ -459,8 +460,10 @@ dft_flux fields::add_dft_flux(const volume_list *where_, const double *freq, siz for (int i = 0; i < 2; ++i) { E = add_dft(cE[i], where->v, freq, Nfreq, true, - where->weight * double(1 - 2 * i), E, false, std::complex(1.0,0), centered_grid); - H = add_dft(cH[i], where->v, freq, Nfreq, false, 1.0, H, false, std::complex(1.0,0), centered_grid); + where->weight * double(1 - 2 * i), E, false, std::complex(1.0,0), + centered_grid, 0, decimation_factor); + H = add_dft(cH[i], where->v, freq, Nfreq, false, 1.0, H, false, std::complex(1.0,0), + centered_grid, 0, decimation_factor); } where = where->next; @@ -544,7 +547,8 @@ double *dft_energy::total() { return F; } -dft_energy fields::add_dft_energy(const volume_list *where_, const double *freq, size_t Nfreq) { +dft_energy fields::add_dft_energy(const volume_list *where_, const double *freq, size_t Nfreq, + int decimation_factor) { if (!where_) // handle empty list of volumes return dft_energy(NULL, NULL, NULL, NULL, freq, Nfreq, v); @@ -555,10 +559,14 @@ dft_energy fields::add_dft_energy(const volume_list *where_, const double *freq, volume_list *where_save = where; while (where) { LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { - E = add_dft(direction_component(Ex, d), where->v, freq, Nfreq, true, 1.0, E); - D = add_dft(direction_component(Dx, d), where->v, freq, Nfreq, false, 1.0, D); - H = add_dft(direction_component(Hx, d), where->v, freq, Nfreq, true, 1.0, H); - B = add_dft(direction_component(Bx, d), where->v, freq, Nfreq, false, 1.0, B); + E = add_dft(direction_component(Ex, d), where->v, freq, Nfreq, true, 1.0, E, + false, 1.0, true, 0, decimation_factor); + D = add_dft(direction_component(Dx, d), where->v, freq, Nfreq, false, 1.0, D, + false, 1.0, true, 0, decimation_factor); + H = add_dft(direction_component(Hx, d), where->v, freq, Nfreq, true, 1.0, H, + false, 1.0, true, 0, decimation_factor); + B = add_dft(direction_component(Bx, d), where->v, freq, Nfreq, false, 1.0, B, + false, 1.0, true, 0, decimation_factor); } where = where->next; } @@ -649,17 +657,18 @@ direction fields::normal_direction(const volume &where) const { } dft_flux fields::add_dft_flux(direction d, const volume &where, const double *freq, size_t Nfreq, - bool use_symmetry, bool centered_grid) { + bool use_symmetry, bool centered_grid, int decimation_factor) { if (d == NO_DIRECTION) d = normal_direction(where); volume_list vl(where, direction_component(Sx, d)); - dft_flux flux = add_dft_flux(&vl, freq, Nfreq, use_symmetry, centered_grid); + dft_flux flux = add_dft_flux(&vl, freq, Nfreq, use_symmetry, centered_grid, decimation_factor); flux.normal_direction = d; return flux; } -dft_flux fields::add_mode_monitor(direction d, const volume &where, const double *freq, size_t Nfreq, bool centered_grid) { - return add_dft_flux(d, where, freq, Nfreq, /*use_symmetry=*/false, centered_grid); +dft_flux fields::add_mode_monitor(direction d, const volume &where, const double *freq, + size_t Nfreq, bool centered_grid, int decimation_factor) { + return add_dft_flux(d, where, freq, Nfreq, /*use_symmetry=*/false, centered_grid, decimation_factor); } dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, diff --git a/src/meep.hpp b/src/meep.hpp index 1479ecfc4..7d58c446c 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1960,25 +1960,35 @@ class fields { bool include_dV = true); 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); + bool use_symmetry = true, bool centered_grid = true, + int decimation_factor = 1); dft_flux add_dft_flux(const volume_list *where, const std::vector &freq, - bool use_symmetry = true, bool centered_grid = true) { - return add_dft_flux(where, freq.data(), freq.size(), use_symmetry, centered_grid); + int decimation_factor = 1, 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) { - return add_dft_flux(where, linspace(freq_min, freq_max, Nfreq), use_symmetry, centered_grid); + bool use_symmetry = true, bool centered_grid = true, + int decimation_factor = 1) { + 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) { - return add_dft_flux(d, where, linspace(freq_min, freq_max, Nfreq), use_symmetry, centered_grid); + int Nfreq, bool use_symmetry = true, bool centered_grid = true, + int decimation_factor = 1) { + 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) { - return add_dft_flux(d, where, freq.data(), freq.size(), use_symmetry, centered_grid); + bool use_symmetry = true, bool centered_grid = true, + int decimation_factor = 1) { + 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); + bool use_symmetry = true, bool centered_grid = true, + int decimation_factor = 1); 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); @@ -1986,13 +1996,17 @@ 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) { - return add_mode_monitor(d, where, linspace(freq_min, freq_max, Nfreq), centered_grid); + int Nfreq, bool centered_grid = true, int decimation_factor = 1) { + 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) { - return add_mode_monitor(d, where, freq.data(), freq.size(), centered_grid); + dft_flux add_mode_monitor(direction d, const volume &where, const std::vector &freq, + bool centered_grid = true, int decimation_factor = 1) { + 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); + dft_flux add_mode_monitor(direction d, const volume &where, const double *freq, size_t Nfreq, + bool centered_grid = true, int decimation_factor = 1); dft_fields add_dft_fields(component *components, int num_components, const volume where, double freq_min, double freq_max, int Nfreq, @@ -2053,34 +2067,42 @@ class fields { void get_mode_mode_overlap(void *mode1_data, void *mode2_data, dft_flux flux, std::complex overlaps[2]); - dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq) { - return add_dft_energy(where, linspace(freq_min, freq_max, Nfreq)); + dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq, + int decimation_factor = 1) { + 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) { - return add_dft_energy(where, freq.data(), freq.size()); + dft_energy add_dft_energy(const volume_list *where, const std::vector &freq, + int decimation_factor = 1) { + 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); + dft_energy add_dft_energy(const volume_list *where, const double *freq, size_t Nfreq, + int decimation_factor = 1); // stress.cpp - dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq) { - return add_dft_force(where, linspace(freq_min, freq_max, Nfreq)); + dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq, + int decimation_factor = 1) { + 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) { - return add_dft_force(where, freq.data(), freq.size()); + dft_force add_dft_force(const volume_list *where, const std::vector &freq, + int decimation_factor = 1) { + 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); + dft_force add_dft_force(const volume_list *where, const double *freq, size_t Nfreq, + int decimation_factor = 1); // near2far.cpp dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, - int Nfreq, int Nperiods = 1) { - return add_dft_near2far(where, linspace(freq_min, freq_max, Nfreq), Nperiods); + int Nfreq, int decimation_factor = 1, 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 Nperiods = 1) { - return add_dft_near2far(where, freq.data(), freq.size(), Nperiods); + int decimation_factor = 1, 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 Nperiods = 1); + int decimation_factor = 1, 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/near2far.cpp b/src/near2far.cpp index 4f097772f..943584bf0 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -557,7 +557,8 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution) static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fabs(a) + fabs(b)); } dft_near2far fields::add_dft_near2far(const volume_list *where, const double *freq, size_t Nfreq, - int Nperiods) { + int decimation_factor, int Nperiods) { + dft_chunk *F = 0; /* E and H chunks*/ double eps = 0, mu = 0; volume everywhere = where->v; @@ -631,7 +632,7 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, const double *fr double s = j == 0 ? 1 : -1; /* sign of n x c */ if (is_electric(c)) s = -s; - F = add_dft(c, w->v, freq, Nfreq, true, s * w->weight, F, false, 1.0, false, c0); + F = add_dft(c, w->v, freq, Nfreq, true, s * w->weight, F, false, 1.0, false, c0, decimation_factor); } } } diff --git a/src/stress.cpp b/src/stress.cpp index 21b4cda36..8f1623d2d 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -151,7 +151,8 @@ void dft_force::scale_dfts(complex scale) { /* note that the components where->c indicate the direction of the force to be computed, so they should be vector components (such as Ex, Ey, ... or Sx, ...) rather than pseudovectors (like Hx, ...). */ -dft_force fields::add_dft_force(const volume_list *where_, const double *freq, size_t Nfreq) { +dft_force fields::add_dft_force(const volume_list *where_, const double *freq, size_t Nfreq, + int decimation_factor) { dft_chunk *offdiag1 = 0, *offdiag2 = 0, *diag = 0; volume_list *where = S.reduce(where_); @@ -167,19 +168,21 @@ dft_force fields::add_dft_force(const volume_list *where_, const double *freq, s if (fd != nd) { // off-diagaonal stress-tensor terms offdiag1 = add_dft(direction_component(Ex, fd), where->v, freq, Nfreq, true, where->weight, - offdiag1); - offdiag2 = add_dft(direction_component(Ex, nd), where->v, freq, Nfreq, false, 1.0, offdiag2); + offdiag1, false, 1.0, true, 0, decimation_factor); + offdiag2 = add_dft(direction_component(Ex, nd), where->v, freq, Nfreq, false, 1.0, offdiag2, + false, 1.0, true, 0, decimation_factor); offdiag1 = add_dft(direction_component(Hx, fd), where->v, freq, Nfreq, true, where->weight, - offdiag1); - offdiag2 = add_dft(direction_component(Hx, nd), where->v, freq, Nfreq, false, 1.0, offdiag2); + offdiag1, false, 1.0, true, 0, decimation_factor); + offdiag2 = add_dft(direction_component(Hx, nd), where->v, freq, Nfreq, false, 1.0, offdiag2, + false, 1.0, true, 0, decimation_factor); } else // diagonal stress-tensor terms LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { complex weight1 = where->weight * (d == fd ? +0.5 : -0.5); diag = add_dft(direction_component(Ex, d), where->v, freq, Nfreq, true, 1.0, diag, true, - weight1, false); + weight1, false, 0, decimation_factor); diag = add_dft(direction_component(Hx, d), where->v, freq, Nfreq, true, 1.0, diag, true, - weight1, false); + weight1, false, 0, decimation_factor); } everywhere = everywhere | where->v; }