diff --git a/doc/docs/Build_From_Source.md b/doc/docs/Build_From_Source.md index f81660d9f..a687b22ad 100644 --- a/doc/docs/Build_From_Source.md +++ b/doc/docs/Build_From_Source.md @@ -257,7 +257,7 @@ make ### Building From Source -The following instructions are for building parallel PyMeep with all optional features from source on Ubuntu 16.04. The parallel version can still be run serially by running a script with just `python` instead of `mpirun -np 4 python`. If you really don't want to install MPI and parallel HDF5, just replace `libhdf5-openmpi-dev` with `libhdf5-dev`, and remove the `--with-mpi`, `CC=mpicc`, and `CPP=mpicxx` flags. The paths to HDF5 will also need to be adjusted to `/usr/lib/x86_64-linux-gnu/hdf5/serial` and `/usr/include/hdf5/serial`. Note that this script builds with Python 3 by default. If you want to use Python 2, just point the `PYTHON` variable to the appropriate interpreter when calling `autogen.sh` for building Meep, and use `pip` instead of `pip3`. +The following instructions are for building parallel PyMeep with all optional features from source on Ubuntu 16.04. (There is a separate [script](http://ab-initio.mit.edu/~oskooi/meep_discuss/build_meep.sh) if you only want the Scheme interface.) The parallel version can still be run serially by running a script with just `python` instead of `mpirun -np 4 python`. If you really don't want to install MPI and parallel HDF5, just replace `libhdf5-openmpi-dev` with `libhdf5-dev`, and remove the `--with-mpi`, `CC=mpicc`, and `CPP=mpicxx` flags. The paths to HDF5 will also need to be adjusted to `/usr/lib/x86_64-linux-gnu/hdf5/serial` and `/usr/include/hdf5/serial`. Note that this script builds with Python 3 by default. If you want to use Python 2, just point the `PYTHON` variable to the appropriate interpreter when calling `autogen.sh` for building Meep, and use `pip` instead of `pip3`. The entire build and install procedure can also be performed using an automated script: diff --git a/doc/docs/Download.md b/doc/docs/Download.md index 8df84d3a5..bd1d72c3b 100644 --- a/doc/docs/Download.md +++ b/doc/docs/Download.md @@ -7,9 +7,7 @@ GitHub Source Repository ------------------------ -The source repository is available via [GitHub](https://github.com/NanoComp/meep). - -The current stable release can be obtained from: +The [source repository](https://github.com/NanoComp/meep) is hosted on GitHub along with gzipped tarballs of stable releases: - @@ -22,7 +20,7 @@ To receive notifications when new versions are released, subscribe to the **meep Precompiled Packages for Ubuntu ------------------------------- -Precompiled packages of Meep version 1.3 (September 2017) are available for [Ubuntu](https://packages.ubuntu.com/search?keywords=meep). We recommend Ubuntu as Meep and all of its dependencies can be installed using just one line: +Precompiled packages of Meep [version 1.3](https://github.com/NanoComp/meep/releases/tag/1.3) (September 2017) are available for [Ubuntu](https://packages.ubuntu.com/search?keywords=meep). We recommend Ubuntu as Meep and all of its dependencies can be installed using just one line: ```sh sudo apt-get install meep h5utils diff --git a/doc/docs/FAQ.md b/doc/docs/FAQ.md index b9c9d8692..6308b4a05 100644 --- a/doc/docs/FAQ.md +++ b/doc/docs/FAQ.md @@ -15,11 +15,11 @@ Meep is a [free and open-source](https://en.wikipedia.org/wiki/Free_and_open-sou ### Who are the developers of Meep? -Meep was originally developed as part of graduate research at MIT. The project is now being maintained by [Simpetus](http://www.simpetus.com) and the developer community on [GitHub](https://github.com/NanoComp/meep). +Meep was originally developed as part of graduate research at MIT. The project has been under continuous development for nearly 20 years. It is currently maintained by [Simpetus](http://www.simpetus.com) and the developer community on [GitHub](https://github.com/NanoComp/meep). ### Where can I ask questions regarding Meep? -There is a public [mailing list](http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss) for users to discuss issues pertaining to setting up simulations, post-processing output, installation, etc. A useful place to start is the [list archives](https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/) which includes all postings (6500+) since 2006 spanning a variety of topics. Bug reports and new feature requests should be filed as a [GitHub issue](https://github.com/NanoComp/meep/issues). +There is a public [mailing list](http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss) for users to discuss issues pertaining to setting up simulations, post-processing output, installation, etc. A useful place to start is the [list archives](https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/) which includes all postings (6500+) since 2006 spanning a variety of topics. Bug reports and new feature requests should be filed as a [GitHub issue](https://github.com/NanoComp/meep/issues). However, do not use issues as a general help desk if you do not understand something (use the mailing list instead). ### How can I contribute to the Meep project? @@ -31,7 +31,7 @@ Yes. The technical details of Meep's inner workings are described in the peer-re ### Where can I find a list of projects which have used Meep? -For a list of more than 2500 published works which have used Meep, see the [Google Scholar citation page](https://scholar.google.com/scholar?hl=en&q=meep+software) as well as that for the [technical reference](https://scholar.google.com/scholar?cites=17712807607104508775) and also the [subpixel smoothing reference](https://scholar.google.com/scholar?cites=410731148689673259). +For a list of more than 2500 published works which have used Meep, see the [Google Scholar citation page](https://scholar.google.com/scholar?hl=en&q=meep+software) as well as that for the [Meep manuscript](https://scholar.google.com/scholar?cites=17712807607104508775) and the [subpixel smoothing reference](https://scholar.google.com/scholar?cites=410731148689673259). For examples based on technology applications, see [simpetus.com/projects.html](http://www.simpetus.com/projects.html). ### Can I access Meep in the public cloud? @@ -54,7 +54,7 @@ Yes. For Windows 10, you can install the [Ubuntu terminal](https://www.microsoft ### Are there precompiled binary packages for Ubuntu? -Yes. Ubuntu and Debian packages can be obtained via the package manager [APT](https://en.wikipedia.org/wiki/APT_(Debian)) as described in [Download](Download.md#precompiled-packages-for-ubuntu). However, the Meep packages for Ubuntu 16.04 ([serial](https://packages.ubuntu.com/xenial/meep) and [parallel](https://packages.ubuntu.com/xenial/meep-openmpi)) and 18.04 ([serial](https://packages.ubuntu.com/bionic/meep) and [parallel](https://packages.ubuntu.com/bionic/meep-openmpi)) are for [version 1.3](https://github.com/NanoComp/meep/releases) (September 2017) which is out of date. The Meep package for Ubuntu is in the process of being updated and will likely appear in Ubuntu 19.10 as derived from the [unstable Debian package](https://packages.debian.org/unstable/meep). In the meantime, since the [Scheme interface](Scheme_User_Interface.md) is no longer being supported and has been replaced by the [Python interface](Python_User_Interface.md), you can use the [Conda packages](Installation.md#conda-packages) which contain the official releases as well as nightly builds of the master branch of the source repository. +Yes. Ubuntu and Debian packages can be obtained via the package manager [APT](https://en.wikipedia.org/wiki/APT_(Debian)) as described in [Download](Download.md#precompiled-packages-for-ubuntu). However, the Meep packages for Ubuntu 16.04 ([serial](https://packages.ubuntu.com/xenial/meep) and [parallel](https://packages.ubuntu.com/xenial/meep-openmpi)) and 18.04 ([serial](https://packages.ubuntu.com/bionic/meep) and [parallel](https://packages.ubuntu.com/bionic/meep-openmpi)) are for [version 1.3](https://github.com/NanoComp/meep/releases/tag/1.3) (September 2017) which is out of date. The Meep package for Ubuntu is in the process of being updated and will likely appear in Ubuntu 19.10 as derived from the [unstable Debian package](https://packages.debian.org/unstable/meep). In the meantime, since the [Scheme interface](Scheme_User_Interface.md) is no longer being supported and has been replaced by the [Python interface](Python_User_Interface.md), you can use the [Conda packages](Installation.md#conda-packages) which contain the official releases as well as nightly builds of the master branch of the source repository. ### Guile is installed, but configure complains that it can't find `guile` @@ -180,7 +180,7 @@ The field from a point source is singular — it blows up as you approach th ### Is a narrow-bandwidth Gaussian pulse considered the same as a continuous-wave (CW) source? -No. A narrow-bandwidth Gaussian is still a Gaussian: it goes to zero at both the beginning and end of its time profile unlike a continuous-wave (CW) source which oscillates indefinitely (but has a [finite turn-on](#why-doesnt-the-continuous-wave-cw-source-produce-an-exact-single-frequency-response)). Assuming you have linear materials, you should get the same results if you use a narrow- or broad-band pulse and look at a single frequency component of the Fourier transform via e.g. [`dft_fields`](Python_User_Interface.md#field-computations). The latter has the advantage that it requires a shorter simulation for the fields to die away due to the [Fourier Uncertainty Principle](https://en.wikipedia.org/wiki/Fourier_transform#Uncertainty_principle). Note also that an almost *zero*-bandwidth Gaussian will produce high-frequency spectral components due to its abrupt turn on and off which are poorly absorbed by PML. +No. A narrow-bandwidth Gaussian is still a Gaussian: it goes to zero at both the beginning and end of its time profile unlike a continuous-wave (CW) source which oscillates indefinitely (but has a [finite turn-on](#why-doesnt-the-continuous-wave-cw-source-produce-an-exact-single-frequency-response)). Assuming you have linear materials, you should get the same results if you use a narrow- or broad-band pulse and look at a single frequency component of the Fourier transform via e.g. [`dft_fields`](Python_User_Interface.md#field-computations). The latter has the advantage that it requires a shorter simulation for the fields to decay away due to the [Fourier Uncertainty Principle](https://en.wikipedia.org/wiki/Fourier_transform#Uncertainty_principle). Note also that an almost *zero*-bandwidth Gaussian will produce high-frequency spectral components due to its abrupt turn on and off which are poorly absorbed by PML. ### How do I create a chirped pulse? @@ -424,7 +424,7 @@ Usage: Other ### Is there a Python interface? -Yes. An official [Python interface](Python_User_Interface.md) was released in [version 1.4](https://github.com/NanoComp/meep/releases) and replaces the [Scheme interface](Scheme_User_Interface.md) which is no longer being supported. An unofficial [Python interface](https://www.fzu.cz/~dominecf/meep/), which predates and is **incompatible** with the official version, has been developed independently by researchers at the Institute of Physics at the Czech Academy of Sciences and Ghent University, and maintained by [Filip Dominec](https://github.com/FilipDominec/python-meep-utils). Unfortunately, this interface has several shortcomings including missing support for geometric objects, lack of high-level abstractions for low-level functionality, and limited documentation. The official interface addresses all these issues. +Yes. An official [Python interface](Python_User_Interface.md) was released in [version 1.4](https://github.com/NanoComp/meep/releases/tag/v1.4.0) and replaces the [Scheme interface](Scheme_User_Interface.md) which is no longer being supported. An unofficial [Python interface](https://www.fzu.cz/~dominecf/meep/), which predates and is **incompatible** with the official version, has been developed independently by researchers at the Institute of Physics at the Czech Academy of Sciences and Ghent University, and maintained by [Filip Dominec](https://github.com/FilipDominec/python-meep-utils). Unfortunately, this interface has several shortcomings including missing support for geometric objects, lack of high-level abstractions for low-level functionality, and limited documentation. The official interface addresses all these issues. ### What is a good rule of thumb for the grid resolution? diff --git a/doc/docs/Field_Functions.md b/doc/docs/Field_Functions.md index b853d5938..3e9900c91 100644 --- a/doc/docs/Field_Functions.md +++ b/doc/docs/Field_Functions.md @@ -128,11 +128,26 @@ meep.Simulation.run(meep.synchronized_magnetic(meep.at_every(1,my_weird_output)) (run-until 200 (synchronized-magnetic (at-every 1 my-weird-output))) ``` -As a final example, the Python interface routine [`get_array_metadata`](Python_User_Interface.md#array-metadata) used to obtain the coordinates of grid points in a volume slice can be replicated in Scheme via e.g.: +As a final example, the function `integrate_field_function` (Python) or `integrate-field-function` (Scheme) can be used to obtain the coordinates of the Yee grid. As long as the field arguments are on the *same* grid (e.g., $E_x$ and $D_x$, $E_y$ and $D_y$, etc.), the integral is computed over the exact Yee grid coordinates rather than being interpolated to the center of each grid point if fields from *different* grids are used (consistent with Meep's paradigm of [pervasive interpolation](Introduction.md#the-illusion-of-continuity)). +**Python** +```py +def f(r,fc): + print("({:.5f}, {:.5f}, {:.5f})".format(r.x,r.y,r.z)) + return 0 + +meep.Simulation.integrate_field_function([mp.Ex],f) +``` + +**Scheme** ```scm -(define (f r eps) (vector3-x r)) -(output-real-field-function "x" Dielectric f) +(define (f r fc) + (begin + (print "(" (vector3-x r) ", " (vector3-y r) ", " (vector3-z r) ")\n") + 0)) +(integrate-field-function (list Ex) f) ``` +This function prints the $(x,y,z)$ Yee grid coordinates of all $E_x$ fields and returns a value of 0 which is never used. In contrast, the output functions `output_field_function` (Python) or `output-field-function` (Scheme) (as well as `output-real-field-function`) interpolate *all* fields onto the center of each grid point. + For more information, see [Python User Interface/Writing Your Own Step Functions](Python_User_Interface.md#writing-your-own-step-functions) or [Scheme User Interface/Writing Your Own Step Functions](Scheme_User_Interface.md#writing-your-own-step-functions). diff --git a/doc/docs/Installation.md b/doc/docs/Installation.md index 19d744e26..0585ce13c 100644 --- a/doc/docs/Installation.md +++ b/doc/docs/Installation.md @@ -22,7 +22,7 @@ Conda Packages The **recommended** way to install PyMeep is using the [Conda](https://conda.io/docs/) package manager. The precompiled binaries run as *fast or faster* than the typical build from source, are simple to install, can be upgraded easily, and take advantage of newer compilers and dependencies than those available in typical systems (e.g., [gcc](https://en.wikipedia.org/wiki/GNU_Compiler_Collection) 7.3.0 vs. Ubuntu 16.04's 5.4; a gap of nearly three years of advances in compiler optimization). Obviously, building from source can still provide advantages if you have access to special hardware or performance libraries that require specific compiler flags (e.g., [icc](https://en.wikipedia.org/wiki/Intel_C%2B%2B_Compiler)); building from source is also required if you are interested in working on the Meep [source code](https://github.com/NanoComp/meep), are performing system-wide installations on a server, or are using systems unsupported by Conda (e.g., supercomputers with Cray MPI). -Binary packages for serial and parallel PyMeep on Linux and macOS are currently available (64 bit architectures only), and are [updated with each new Meep release](https://github.com/conda-forge/pymeep-feedstock). (Note: the Conda packages will *not* work on [Windows](#installation-on-windows).) The easiest way to get started is to install [Miniconda](https://conda.io/miniconda.html), which comes with everything necessary to create Python environments with Conda. For example, to install Miniconda with Python 3 on Linux: +Binary packages for serial and parallel PyMeep on Linux and macOS are currently available (64 bit architectures only), and are [updated with each new Meep release](https://github.com/conda-forge/pymeep-feedstock). Note: the Conda packages will *not* work on native [Windows](#installation-on-windows) (unless you install the Ubuntu terminal app) and do *not* include the Scheme interface which must be [built from source](Build_From_Source.md). The easiest way to get started is to install [Miniconda](https://conda.io/miniconda.html), which comes with everything necessary to create Python environments with Conda. For example, to install Miniconda with Python 3 on Linux: ```bash wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh @@ -171,4 +171,4 @@ To build the latest version of Meep from source on macOS Sierra, follow these [i Installation on Windows ---------------------------- -Native Windows installation is currently unsupported. The recommended procedure is to install Ubuntu using the [Windows Subsystem for Linux](https://docs.microsoft.com/en-us/windows/wsl/install-win10). This gives you access to a bash terminal running on Ubuntu from within Windows. From there you can install the Conda packages as described above. The drawback is that you can't see plots from matplotlib (though saving them to disk and opening them from Windows works fine). The easiest way around this is to add the `jupyter` package to the `conda create ...` command. This will allow you to run a [Jupyter notebook](https://jupyter.readthedocs.io/en/latest/) in the browser, and from there you can visualize plots interactively. +Native Windows installation is currently unsupported. The recommended procedure is to install Ubuntu using the [Windows Subsystem for Linux](https://docs.microsoft.com/en-us/windows/wsl/install-win10). This gives you access to a bash terminal running Ubuntu from within Windows. From there you can install the Conda packages as described above. The drawback is that you can't see plots from matplotlib (though saving them to disk and opening them from Windows works fine). The easiest way around this is to add the `jupyter` package to the `conda create ...` command. This will allow you to run a [Jupyter notebook](https://jupyter.readthedocs.io/en/latest/) in the browser, and from there you can visualize plots interactively. diff --git a/doc/docs/Introduction.md b/doc/docs/Introduction.md index a3bcaa5fc..be8f1957c 100644 --- a/doc/docs/Introduction.md +++ b/doc/docs/Introduction.md @@ -6,7 +6,7 @@ Meep implements the [finite-difference time-domain](https://en.wikipedia.org/wik This section introduces the equations and the electromagnetic units employed by Meep, the FDTD method, and Meep's approach to FDTD. Also, FDTD is only one of several useful methods in computational electromagnetics, each of which has their own special uses — a few of the other methods are mentioned, and some hints are provided as to which applications FDTD is well suited for and when you should potentially consider a different method. -This introduction does not describe the [User Interface](Python_User_Interface.md) with which you set up simulations. Instead, the focus here is on the physics and numerical methods. For tutorial examples which demonstrate core functionality, see [Tutorial/Basics](Python_Tutorials/Basics.md). +This introduction does not describe the [Python User Interface](Python_User_Interface.md) with which you set up simulations. Instead, the focus here is on the physics and numerical methods. For tutorial examples which demonstrate core functionality, see [Tutorial/Basics](Python_Tutorials/Basics.md) and the [Simpetus projects page](http://www.simpetus.com/projects.html). [TOC] @@ -100,7 +100,7 @@ FDTD is, of course, not the only numerical method in computational electromagnet For example, although FDTD can be used to compute electromagnetic eigenmodes (below), in lossless structures it is often quicker, easier, and more reliable to use a specialized eigenmode solver such as [MPB](http://mpb.readthedocs.io). See also the [frequency vs. time domain](http://mpb.readthedocs.io/en/latest/Introduction/) discussion in the MPB manual and the [resonant modes](#resonant-modes) discussion below. -For computing the field pattern or response of a structure at a *single frequency*, it may be more efficient to directly solve the corresponding linear equation rather than iterating in time. Indeed, this can be done directly in Meep (i.e. a [finite-difference frequency-domain solver](Python_User_Interface.md#frequency-domain-solver)) — see [Tutorials/Frequency-Domain Solver](Python_Tutorials/Frequency_Domain_Solver.md). However, especially in cases where there are large differences in scale (e.g. with metals with a shallow skin depth), it may be better to use a method that allows a variable resolution in different spatial regions, such as a finite-element or boundary-element method. Boundary-element methods are especially powerful when you have a large volume-to-surface ratio, such as for scattering calculations over small objects in a large (i.e., infinite-sized) volume. +For computing the field pattern or response of a structure at a *single frequency*, it may be more efficient to directly solve the corresponding linear equation rather than iterating in time. Indeed, this can be done directly in Meep (i.e. a [finite-difference frequency-domain solver](Python_User_Interface.md#frequency-domain-solver)) — see [Tutorial/Frequency-Domain Solver](Python_Tutorials/Frequency_Domain_Solver.md). However, especially in cases where there are large differences in scale (e.g. with metals with a shallow skin depth), it may be better to use a method that allows a variable resolution in different spatial regions, such as a finite-element or boundary-element method. Boundary-element methods are especially powerful when you have a large volume-to-surface ratio, such as for scattering calculations over small objects in a large (i.e., infinite-sized) volume. A strength of time-domain methods is their ability to obtain the [entire frequency spectrum of responses (or eigenfrequencies) in a single simulation](Python_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend), by Fourier-transforming the response to a short pulse or using more sophisticated signal-processing methods such as [Harminv](Python_User_Interface.md#harminv). Finite-element methods can also be used for time-evolving fields, but they suffer a serious disadvantage compared to finite-difference methods: finite-element methods, for stability, must typically use some form of *implicit time-stepping*, where they must invert a matrix (solve a linear system) at every time step. diff --git a/doc/docs/index.md b/doc/docs/index.md index d7128919f..c25317b32 100644 --- a/doc/docs/index.md +++ b/doc/docs/index.md @@ -41,7 +41,7 @@ Meep's scriptable interface makes it possible to combine many sorts of computati Download -------- -The source repository is on [GitHub](https://github.com/NanoComp/meep). Gzipped tarballs of stable versions are in [Releases](https://github.com/NanoComp/meep/releases). The release history is in [NEWS](https://github.com/NanoComp/meep/blob/master/NEWS.md). Installation instructions are in [Installation](Installation.md). +The [source repository](https://github.com/NanoComp/meep) is hosted on GitHub. Gzipped tarballs of stable versions are in [Releases](https://github.com/NanoComp/meep/releases). The release history is in [NEWS](https://github.com/NanoComp/meep/blob/master/NEWS.md). Installation instructions are in [Installation](Installation.md). Documentation ------------- diff --git a/python/meep.i b/python/meep.i index 06e4ad345..0024791cd 100644 --- a/python/meep.i +++ b/python/meep.i @@ -307,7 +307,7 @@ PyObject *py_do_harminv(PyObject *vals, double dt, double f_min, double f_max, i // Wrapper around meep::dft_near2far::farfield PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) { - Py_ssize_t len = f->Nfreq * 6; + Py_ssize_t len = f->freq.size() * 6; PyObject *res = PyList_New(len); std::complex *ff_arr = f->farfield(v); @@ -336,7 +336,7 @@ PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) { if (!EH) return PyArray_SimpleNew(0, 0, NPY_CDOUBLE); // frequencies are the last dimension - if (n2f->Nfreq > 1) dims[rank++] = n2f->Nfreq; + if (n2f->freq.size() > 1) dims[rank++] = n2f->freq.size(); // Additional rank to store all 12 E/H x/y/z r/i arrays. rank++; @@ -347,7 +347,7 @@ PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) { } PyObject *py_arr = PyArray_SimpleNew(rank, arr_dims, NPY_DOUBLE); - memcpy(PyArray_DATA((PyArrayObject*)py_arr), EH, sizeof(meep::realnum) * 2 * N * 6 * n2f->Nfreq); + memcpy(PyArray_DATA((PyArrayObject*)py_arr), EH, sizeof(meep::realnum) * 2 * N * 6 * n2f->freq.size()); delete[] arr_dims; delete[] EH; @@ -357,7 +357,7 @@ PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) { // Wrapper around meep::dft_ldos::ldos PyObject *_dft_ldos_ldos(meep::dft_ldos *f) { - Py_ssize_t len = f->Nomega; + Py_ssize_t len = f->freq.size(); PyObject *res = PyList_New(len); double *tmp = f->ldos(); @@ -373,7 +373,7 @@ PyObject *_dft_ldos_ldos(meep::dft_ldos *f) { // Wrapper around meep::dft_ldos_F PyObject *_dft_ldos_F(meep::dft_ldos *f) { - Py_ssize_t len = f->Nomega; + Py_ssize_t len = f->freq.size(); PyObject *res = PyList_New(len); std::complex *tmp = f->F(); @@ -389,7 +389,7 @@ PyObject *_dft_ldos_F(meep::dft_ldos *f) { // Wrapper arond meep::dft_ldos_J PyObject *_dft_ldos_J(meep::dft_ldos *f) { - Py_ssize_t len = f->Nomega; + Py_ssize_t len = f->freq.size(); PyObject *res = PyList_New(len); std::complex *tmp = f->J(); @@ -455,7 +455,7 @@ void _get_dft_data(meep::dft_chunk *dc, std::complex *cdata, int } for (meep::dft_chunk *cur = dc; cur; cur = cur->next_in_dft) { - size_t Nchunk = cur->N * cur->Nomega; + size_t Nchunk = cur->N * cur->omega.size(); for (size_t i = 0; i < Nchunk; ++i) { cdata[i + istart] = cur->dft[i]; } @@ -473,7 +473,7 @@ void _load_dft_data(meep::dft_chunk *dc, std::complex *cdata, int } for (meep::dft_chunk *cur = dc; cur; cur = cur->next_in_dft) { - size_t Nchunk = cur->N * cur->Nomega; + size_t Nchunk = cur->N * cur->omega.size(); for (size_t i = 0; i < Nchunk; ++i) { cur->dft[i] = cdata[i + istart]; } @@ -494,7 +494,7 @@ kpoint_list get_eigenmode_coefficients_and_kpoints(meep::fields *f, meep::dft_fl double *vgrp, meep::kpoint_func user_kpoint_func, void *user_kpoint_data, double *cscale, meep::direction d) { - size_t num_kpoints = num_bands * flux.Nfreq; + size_t num_kpoints = num_bands * flux.freq.size(); meep::vec *kpoints = new meep::vec[num_kpoints]; meep::vec *kdom = new meep::vec[num_kpoints]; @@ -800,7 +800,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, // Typemap suite for dft_flux %typemap(out) double* flux { - int size = arg1->Nfreq; + int size = arg1->freq.size(); $result = PyList_New(size); for(int i = 0; i < size; i++) { PyList_SetItem($result, i, PyFloat_FromDouble($1[i])); @@ -812,7 +812,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, // Typemap suite for dft_force %typemap(out) double* force { - int size = arg1->Nfreq; + int size = arg1->freq.size(); $result = PyList_New(size); for(int i = 0; i < size; i++) { PyList_SetItem($result, i, PyFloat_FromDouble($1[i])); @@ -1343,35 +1343,6 @@ void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const double kdom[3]); #endif // HAVE_MPB -// Make omega members of meep::dft_ldos available as 'freq' in python -%extend meep::dft_ldos { - - double get_omega_min() { - return $self->omega_min; - } - double get_domega() { - return $self->domega; - } - int get_Nomega() { - return $self->Nomega; - } - %pythoncode %{ - def freqs(self): - import math - import numpy as np - start = self.omega_min / (2 * math.pi) - stop = start + (self.domega / (2 * math.pi)) * self.Nomega - return np.linspace(start, stop, num=self.Nomega, endpoint=False).tolist() - - __swig_getmethods__["freq_min"] = get_omega_min - __swig_getmethods__["nfreq"] = get_Nomega - __swig_getmethods__["dfreq"] = get_domega - if _newclass: freq_min = property(get_omega_min) - if _newclass: nfreq = property(get_Nomega) - if _newclass: dfreq = property(get_domega) - %} -} - %extend meep::fields { bool is_periodic(boundary_side side, direction dir) { return $self->boundaries[side][dir] == meep::Periodic; diff --git a/python/simulation.py b/python/simulation.py index a58e42430..28b95442a 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -316,16 +316,8 @@ def remove(self): return self.swigobj_attr('remove') @property - def freq_min(self): - return self.swigobj_attr('freq_min') - - @property - def dfreq(self): - return self.swigobj_attr('dfreq') - - @property - def Nfreq(self): - return self.swigobj_attr('Nfreq') + def freq(self): + return self.swigobj_attr('freq') @property def where(self): @@ -364,6 +356,9 @@ def cH(self): def normal_direction(self): return self.swigobj_attr('normal_direction') + @property + def freq(self): + return self.swigobj_attr('freq') class DftForce(DftObj): @@ -389,6 +384,9 @@ def offdiag2(self): def diag(self): return self.swigobj_attr('diag') + @property + def freq(self): + return self.swigobj_attr('freq') class DftNear2Far(DftObj): @@ -422,6 +420,9 @@ def mu(self): def flux(self, direction, where, resolution): return self.swigobj_attr('flux')(direction, where.swigobj, resolution) + @property + def freq(self): + return self.swigobj_attr('freq') class DftEnergy(DftObj): @@ -443,6 +444,9 @@ def magnetic(self): def total(self): return self.swigobj_attr('total') + @property + def freq(self): + return self.swigobj_attr('freq') class DftFields(DftObj): @@ -714,15 +718,14 @@ def _get_valid_material_frequencies(self): def _check_material_frequencies(self): min_freq, max_freq = self._get_valid_material_frequencies() - source_freqs = [(s.src.frequency, 0 if s.src.width == 0 else 1 / s.src.width) for s in self.sources if hasattr(s.src, 'frequency')] dft_freqs = [] for dftf in self.dft_objects: - dft_freqs.append(dftf.freq_min) - dft_freqs.append(dftf.freq_min + dftf.Nfreq * dftf.dfreq) + dft_freqs.append(dftf.freq[0]) + dft_freqs.append(dftf.freq[-1]) warn_src = ('Note: your sources include frequencies outside the range of validity of the ' + 'material models. This is fine as long as you eventually only look at outputs ' + @@ -2104,9 +2107,9 @@ def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, eig_v direction = flux.normal_direction num_bands = len(bands) - coeffs = np.zeros(2 * num_bands * flux.Nfreq, dtype=np.complex128) - vgrp = np.zeros(num_bands * flux.Nfreq) - cscale = np.zeros(num_bands * flux.Nfreq) + coeffs = np.zeros(2 * num_bands * flux.freq.size(), dtype=np.complex128) + vgrp = np.zeros(num_bands * flux.freq.size()) + cscale = np.zeros(num_bands * flux.freq.size()) kpoints, kdom = mp.get_eigenmode_coefficients_and_kpoints( self.fields, @@ -2123,7 +2126,7 @@ def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, eig_v direction ) - return EigCoeffsResult(np.reshape(coeffs, (num_bands, flux.Nfreq, 2)), vgrp, kpoints, kdom, cscale) + return EigCoeffsResult(np.reshape(coeffs, (num_bands, flux.freq.size(), 2)), vgrp, kpoints, kdom, cscale) def get_eigenmode(self, freq, direction, where, band_num, kpoint, eig_vol=None, match_frequency=True, parity=mp.NO_PARITY, resolution=0, eigensolver_tol=1e-12): @@ -2958,7 +2961,7 @@ def _ldos(sim, todo): sim.ldos_data = mp._dft_ldos_ldos(ldos) sim.ldos_Fdata = mp._dft_ldos_F(ldos) sim.ldos_Jdata = mp._dft_ldos_J(ldos) - display_csv(sim, 'ldos', zip(ldos.freqs(), sim.ldos_data)) + display_csv(sim, 'ldos', zip(ldos.freq, sim.ldos_data)) return _ldos @@ -2967,7 +2970,7 @@ def scale_flux_fields(s, flux): def get_flux_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return [f.freq[i] for i in range(f.freq.size())] def get_fluxes(f): @@ -2979,11 +2982,11 @@ def scale_force_fields(s, force): def get_eigenmode_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return [f.freq[i] for i in range(f.freq.size())] def get_force_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return [f.freq[i] for i in range(f.freq.size())] def get_forces(f): @@ -2995,7 +2998,7 @@ def scale_near2far_fields(s, n2f): def get_near2far_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return [f.freq[i] for i in range(f.freq.size())] def scale_energy_fields(s, ef): @@ -3003,7 +3006,7 @@ def scale_energy_fields(s, ef): def get_energy_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return [f.freq[i] for i in range(f.freq.size())] def get_electric_energy(f): diff --git a/python/tests/geom.py b/python/tests/geom.py index e23e04ec7..6ac6557e1 100644 --- a/python/tests/geom.py +++ b/python/tests/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(15, 6, 2, fregion) + sim.add_flux(18, 6, 2, fregion) check_warnings(sim) # Invalid geometry material sim = mp.Simulation(cell_size=cell_size, resolution=resolution, sources=valid_sources[0], geometry=geom) - sim.add_flux(15, 6, 2, fregion) + sim.add_flux(18, 6, 2, fregion) check_warnings(sim) def test_transform(self): diff --git a/python/tests/ldos.py b/python/tests/ldos.py index 726e3cd11..481a7fcb5 100644 --- a/python/tests/ldos.py +++ b/python/tests/ldos.py @@ -43,10 +43,7 @@ def test_ldos_user_object(self): ) self.assertAlmostEqual(self.sim.ldos_data[0], 1.011459560620368) - freqs = ldos.freqs() - self.assertEqual(ldos.freq_min, freqs[0] * 2 * math.pi) - self.assertEqual(ldos.nfreq, 1) - self.assertEqual(ldos.dfreq, 0) + self.assertEqual(ldos.freq.size(), 1) def test_invalid_dft_ldos(self): with self.assertRaises(ValueError): diff --git a/scheme/meep.cpp b/scheme/meep.cpp index c8526b236..dbc653e3d 100644 --- a/scheme/meep.cpp +++ b/scheme/meep.cpp @@ -61,7 +61,7 @@ kpoint_list do_get_eigenmode_coefficients(fields *f, dft_flux flux, const volume double eig_resolution, double eigensolver_tol, meep::kpoint_func user_kpoint_func, void *user_kpoint_data, int dir) { - size_t num_kpoints = num_bands * flux.Nfreq; + size_t num_kpoints = num_bands * flux.freq.size(); meep::vec *kpoints = new meep::vec[num_kpoints]; meep::vec *kdom = new meep::vec[num_kpoints]; double *cscale = 0; // Not needed until adjoint calculation is implemented in scheme @@ -91,63 +91,63 @@ volume_list *make_volume_list(const volume &v, int c, complex weight, vo ctlio::number_list dft_flux_flux(dft_flux *f) { ctlio::number_list res; - res.num_items = f->Nfreq; + res.num_items = f->freq.size(); res.items = f->flux(); return res; } ctlio::number_list dft_energy_electric(dft_energy *f) { ctlio::number_list res; - res.num_items = f->Nfreq; + res.num_items = f->freq.size(); res.items = f->electric(); return res; } ctlio::number_list dft_energy_magnetic(dft_energy *f) { ctlio::number_list res; - res.num_items = f->Nfreq; + res.num_items = f->freq.size(); res.items = f->magnetic(); return res; } ctlio::number_list dft_energy_total(dft_energy *f) { ctlio::number_list res; - res.num_items = f->Nfreq; + res.num_items = f->freq.size(); res.items = f->total(); return res; } ctlio::number_list dft_force_force(dft_force *f) { ctlio::number_list res; - res.num_items = f->Nfreq; + res.num_items = f->freq.size(); res.items = f->force(); return res; } ctlio::number_list dft_ldos_ldos(dft_ldos *f) { ctlio::number_list res; - res.num_items = f->Nomega; + res.num_items = f->freq.size(); res.items = f->ldos(); return res; } ctlio::cnumber_list dft_ldos_F(dft_ldos *f) { ctlio::cnumber_list res; - res.num_items = f->Nomega; + res.num_items = f->freq.size(); res.items = (cnumber *)f->F(); return res; } ctlio::cnumber_list dft_ldos_J(dft_ldos *f) { ctlio::cnumber_list res; - res.num_items = f->Nomega; + res.num_items = f->freq.size(); res.items = (cnumber *)f->J(); return res; } ctlio::cnumber_list dft_near2far_farfield(dft_near2far *f, const vec &x) { ctlio::cnumber_list res; - res.num_items = f->Nfreq * 6; + res.num_items = f->freq.size() * 6; res.items = (cnumber *)f->farfield(x); return res; } @@ -155,7 +155,7 @@ ctlio::cnumber_list dft_near2far_farfield(dft_near2far *f, const vec &x) { ctlio::number_list dft_near2far_flux(dft_near2far *f, direction df, const volume &where, double resolution) { ctlio::number_list res; - res.num_items = f->Nfreq; + res.num_items = f->freq.size(); res.items = f->flux(df, where, resolution); return res; } diff --git a/src/dft.cpp b/src/dft.cpp index af5452fd1..28ecaaee7 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -30,11 +30,22 @@ typedef complex cdouble; namespace meep { +std::vector linspace(double freq_min, double freq_max, int Nfreq) { + double dfreq = Nfreq <= 1 ? 0.0 : (freq_max - freq_min) / (Nfreq - 1); + std::vector freq; + if (Nfreq <= 1) + freq.push_back((freq_min + freq_max) * 0.5); + else + for (int i = 0; i < Nfreq; ++i) + freq.push_back(freq_min + i*dfreq); + + return freq; +} + struct dft_chunk_data { // for passing to field::loop_in_chunks as void* component c; int vc; - double omega_min, domega; - int Nomega; + std::vector omega; complex stored_weight, extra_weight; double dt_factor; bool include_dV_and_interp_weights; @@ -84,9 +95,8 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve sn = sn_; vc = data->vc; - omega_min = data->omega_min; - domega = data->domega; - Nomega = data->Nomega; + const int Nomega = data->omega.size(); + omega = data->omega; dft_phase = new complex[Nomega]; N = 1; @@ -150,6 +160,16 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do complex stored_weight, dft_chunk *chunk_next, bool sqrt_dV_and_interp_weights, complex extra_weight, bool use_centered_grid, int vc) { + return add_dft(c, where, meep::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); +} + +dft_chunk *fields::add_dft(component c, const volume &where, const std::vector freq, + 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) { if (coordinate_mismatch(gv.dim, c)) return NULL; /* If you call add_dft before adding sources, it will do nothing @@ -164,10 +184,8 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do dft_chunk_data data; data.c = c; data.vc = vc; - if (Nfreq <= 1) freq_min = freq_max = (freq_min + freq_max) * 0.5; - data.omega_min = freq_min * 2 * pi; - data.domega = Nfreq <= 1 ? 0.0 : (freq_max * 2 * pi - data.omega_min) / (Nfreq - 1); - data.Nomega = Nfreq; + for (int i = 0; i < freq.size(); ++i) + data.omega.push_back(2*pi*freq[i]); data.stored_weight = stored_weight; data.extra_weight = extra_weight; data.dt_factor = dt / sqrt(2.0 * pi); @@ -184,11 +202,16 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool include_dV_and_interp_weights) { + return add_dft(where, meep::linspace(freq_min, freq_max, Nfreq), include_dV_and_interp_weights); +} + +dft_chunk *fields::add_dft(const volume_list *where, const std::vector freq, + bool include_dV_and_interp_weights) { dft_chunk *chunks = 0; while (where) { if (is_derived(where->c)) abort("derived_component invalid for dft"); cdouble stored_weight = where->weight; - chunks = add_dft(component(where->c), where->v, freq_min, freq_max, Nfreq, + chunks = add_dft(component(where->c), where->v, freq, include_dV_and_interp_weights, stored_weight, chunks); where = where->next; } @@ -197,7 +220,11 @@ dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double fre dft_chunk *fields::add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq) { - return add_dft(c, where, freq_min, freq_max, Nfreq, false); + return add_dft(c, where, meep::linspace(freq_min, freq_max, Nfreq), false); +} + +dft_chunk *fields::add_dft_pt(component c, const vec &where, const std::vector freq) { + return add_dft(c, where, freq, false); } void fields::update_dfts() { @@ -217,8 +244,9 @@ void fields_chunk::update_dfts(double timeE, double timeH) { void dft_chunk::update_dft(double time) { if (!fc->f[c][0]) return; + const int Nomega = omega.size(); for (int i = 0; i < Nomega; ++i) - dft_phase[i] = polar(1.0, (omega_min + i * domega) * time) * scale; + dft_phase[i] = polar(1.0, omega[i] * time) * scale; int numcmp = fc->f[c][1] ? 2 : 1; @@ -258,16 +286,16 @@ void dft_chunk::update_dft(double time) { } void dft_chunk::scale_dft(complex scale) { - for (size_t i = 0; i < N * Nomega; ++i) + for (size_t i = 0; i < N * omega.size(); ++i) dft[i] *= scale; if (next_in_dft) next_in_dft->scale_dft(scale); } void dft_chunk::operator-=(const dft_chunk &chunk) { - if (c != chunk.c || N * Nomega != chunk.N * chunk.Nomega) + if (c != chunk.c || N * omega.size() != chunk.N * chunk.omega.size()) abort("Mismatched chunks in dft_chunk::operator-="); - for (size_t i = 0; i < N * Nomega; ++i) + for (size_t i = 0; i < N * omega.size(); ++i) dft[i] -= chunk.dft[i]; if (next_in_dft) { @@ -279,7 +307,7 @@ void dft_chunk::operator-=(const dft_chunk &chunk) { size_t dft_chunks_Ntotal(dft_chunk *dft_chunks, size_t *my_start) { size_t n = 0; for (dft_chunk *cur = dft_chunks; cur; cur = cur->next_in_dft) - n += cur->N * cur->Nomega * 2; + n += cur->N * cur->omega.size() * 2; *my_start = partial_sum_to_all(n) - n; // sum(n) for processes before this return sum_to_all(n); } @@ -297,7 +325,7 @@ void save_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const file->create_data(dataname, 1, &n); for (dft_chunk *cur = dft_chunks; cur; cur = cur->next_in_dft) { - size_t Nchunk = cur->N * cur->Nomega * 2; + size_t Nchunk = cur->N * cur->omega.size() * 2; file->write_chunk(1, &istart, &Nchunk, (realnum *)cur->dft); istart += Nchunk; } @@ -325,7 +353,7 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const file->file_name(), dataname); for (dft_chunk *cur = dft_chunks; cur; cur = cur->next_in_dft) { - size_t Nchunk = cur->N * cur->Nomega * 2; + size_t Nchunk = cur->N * cur->omega.size() * 2; file->read_chunk(1, &istart, &Nchunk, (realnum *)cur->dft); istart += Nchunk; } @@ -338,17 +366,21 @@ void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_) - : Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), + : E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), use_symmetry(use_symmetry_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + freq = meep::linspace(fmin, fmax, Nf); +} + +dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, + const std::vector freq_, const volume &where_, + direction normal_direction_, bool use_symmetry_) + : E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), + use_symmetry(use_symmetry_) { + freq = freq_; } dft_flux::dft_flux(const dft_flux &f) : where(f.where) { - freq_min = f.freq_min; - Nfreq = f.Nfreq; - dfreq = f.dfreq; + freq = f.freq; E = f.E; H = f.H; cE = f.cE; @@ -358,6 +390,7 @@ dft_flux::dft_flux(const dft_flux &f) : where(f.where) { } double *dft_flux::flux() { + const int Nfreq = freq.size(); double *F = new double[Nfreq]; for (int i = 0; i < Nfreq; ++i) F[i] = 0; @@ -403,8 +436,13 @@ void dft_flux::scale_dfts(complex scale) { dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double freq_max, int Nfreq, bool use_symmetry) { + return add_dft_flux(where_, meep::linspace(freq_min, freq_max, Nfreq), use_symmetry); +} + +dft_flux fields::add_dft_flux(const volume_list *where_, const std::vector freq, + bool use_symmetry) { if (!where_) // handle empty list of volumes - return dft_flux(Ex, Hy, NULL, NULL, freq_min, freq_max, Nfreq, v, NO_DIRECTION, use_symmetry); + return dft_flux(Ex, Hy, NULL, NULL, freq, v, NO_DIRECTION, use_symmetry); dft_chunk *E = 0, *H = 0; component cE[2] = {Ex, Ey}, cH[2] = {Hy, Hx}; @@ -437,9 +475,9 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double } for (int i = 0; i < 2; ++i) { - E = add_dft(cE[i], where->v, freq_min, freq_max, Nfreq, true, + E = add_dft(cE[i], where->v, freq, true, where->weight * double(1 - 2 * i), E); - H = add_dft(cH[i], where->v, freq_min, freq_max, Nfreq, false, 1.0, H); + H = add_dft(cH[i], where->v, freq, false, 1.0, H); } where = where->next; @@ -449,21 +487,23 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double // if the volume list has only one entry, store its component's direction. // if the volume list has > 1 entry, store NO_DIRECTION. direction flux_dir = (where_->next ? NO_DIRECTION : component_direction(where_->c)); - return dft_flux(cE[0], cH[0], E, H, freq_min, freq_max, Nfreq, firstvol, flux_dir, use_symmetry); + return dft_flux(cE[0], cH[0], E, H, freq, firstvol, flux_dir, use_symmetry); +} + +dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, + double fmin, double fmax, int Nf, const volume &where_) + : E(E_), H(H_), D(D_), B(B_), where(where_) { + freq = meep::linspace(fmin, fmax, Nf); } -dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, - double fmax, int Nf, const volume &where_) - : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); +dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, + const std::vector freq_, const volume &where_) + : E(E_), H(H_), D(D_), B(B_), where(where_) { + freq = freq_; } dft_energy::dft_energy(const dft_energy &f) : where(f.where) { - freq_min = f.freq_min; - Nfreq = f.Nfreq; - dfreq = f.dfreq; + freq = f.freq; E = f.E; H = f.H; D = f.D; @@ -471,6 +511,7 @@ dft_energy::dft_energy(const dft_energy &f) : where(f.where) { } double *dft_energy::electric() { + const int Nfreq = freq.size(); double *F = new double[Nfreq]; for (int i = 0; i < Nfreq; ++i) F[i] = 0; @@ -486,6 +527,7 @@ double *dft_energy::electric() { } double *dft_energy::magnetic() { + const int Nfreq = freq.size(); double *F = new double[Nfreq]; for (int i = 0; i < Nfreq; ++i) F[i] = 0; @@ -501,6 +543,7 @@ double *dft_energy::magnetic() { } double *dft_energy::total() { + const int Nfreq = freq.size(); double *Fe = electric(); double *Fm = magnetic(); double *F = new double[Nfreq]; @@ -513,9 +556,13 @@ double *dft_energy::total() { dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, int Nfreq) { + return add_dft_energy(where_, meep::linspace(freq_min, freq_max, Nfreq)); +} + +dft_energy fields::add_dft_energy(const volume_list *where_, const std::vector freq) { if (!where_) // handle empty list of volumes - return dft_energy(NULL, NULL, NULL, NULL, freq_min, freq_max, Nfreq, v); + return dft_energy(NULL, NULL, NULL, NULL, freq, v); dft_chunk *E = 0, *D = 0, *H = 0, *B = 0; volume firstvol(where_->v); @@ -523,16 +570,16 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, do volume_list *where_save = where; while (where) { LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { - E = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, E); - D = add_dft(direction_component(Dx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, D); - H = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, H); - B = add_dft(direction_component(Bx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, B); + E = add_dft(direction_component(Ex, d), where->v, freq, true, 1.0, E); + D = add_dft(direction_component(Dx, d), where->v, freq, false, 1.0, D); + H = add_dft(direction_component(Hx, d), where->v, freq, true, 1.0, H); + B = add_dft(direction_component(Bx, d), where->v, freq, false, 1.0, B); } where = where->next; } delete where_save; - return dft_energy(E, H, D, B, freq_min, freq_max, Nfreq, firstvol); + return dft_energy(E, H, D, B, freq, firstvol); } void dft_energy::save_hdf5(h5file *file, const char *dprefix) { @@ -615,20 +662,32 @@ direction fields::normal_direction(const volume &where) const { dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry) { + return add_dft_flux(d, where, meep::linspace(freq_min, freq_max, Nfreq), use_symmetry); +} + +dft_flux fields::add_dft_flux(direction d, const volume &where, const std::vector freq, + bool use_symmetry) { if (d == NO_DIRECTION) d = normal_direction(where); volume_list vl(where, direction_component(Sx, d)); - dft_flux flux = add_dft_flux(&vl, freq_min, freq_max, Nfreq, use_symmetry); + dft_flux flux = add_dft_flux(&vl, freq, use_symmetry); flux.normal_direction = d; return flux; } dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, int Nfreq) { - return add_dft_flux(d, where, freq_min, freq_max, Nfreq, /*use_symmetry=*/false); + return add_mode_monitor(d, where, meep::linspace(freq_min, freq_max, Nfreq)); } -dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, - int Nfreq) { +dft_flux fields::add_mode_monitor(direction d, const volume &where, const std::vector freq) { + return add_dft_flux(d, where, freq, /*use_symmetry=*/false); +} + +dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq) { + return add_dft_flux_box(where, meep::linspace(freq_min, freq_max, Nfreq)); +} + +dft_flux fields::add_dft_flux_box(const volume &where, const std::vector freq) { volume_list *faces = 0; LOOP_OVER_DIRECTIONS(where.dim, d) { if (where.in_direction(d) > 0) { @@ -642,23 +701,29 @@ dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double f } } - dft_flux flux = add_dft_flux(faces, freq_min, freq_max, Nfreq); + dft_flux flux = add_dft_flux(faces, freq); delete faces; return flux; } -dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, - int Nfreq) { - return add_dft_flux(NO_DIRECTION, where, freq_min, freq_max, Nfreq); +dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq) { + return add_dft_flux_plane(where, meep::linspace(freq_min, freq_max, Nfreq)); +} + +dft_flux fields::add_dft_flux_plane(const volume &where, const std::vector freq) { + return add_dft_flux(NO_DIRECTION, where, freq); +} + +dft_fields::dft_fields(dft_chunk *chunks_, double freq_min, double freq_max, int Nf, const volume &where_) + : where(where_) { + chunks = chunks_; + freq = meep::linspace(freq_min, freq_max, Nf); } -dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, int Nfreq_, - const volume &where_) +dft_fields::dft_fields(dft_chunk *chunks_, const std::vector freq_, const volume &where_) : where(where_) { chunks = chunks_; - freq_min = freq_min_; - dfreq = Nfreq_ <= 1 ? 0.0 : (freq_max_ - freq_min_) / (Nfreq_ - 1); - Nfreq = Nfreq_; + freq = freq_; } void dft_fields::scale_dfts(cdouble scale) { chunks->scale_dft(scale); } @@ -673,17 +738,23 @@ void dft_fields::remove() { dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { + return add_dft_fields(components, num_components, where, + meep::linspace(freq_min, freq_max, Nfreq), use_centered_grid); +} + +dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, + const std::vector freq, bool use_centered_grid) { bool include_dV_and_interp_weights = false; bool sqrt_dV_and_interp_weights = false; // default option from meep.hpp (expose to user?) std::complex extra_weight = 1.0; // default option from meep.hpp (expose to user?) cdouble stored_weight = 1.0; dft_chunk *chunks = 0; for (int nc = 0; nc < num_components; nc++) - chunks = add_dft(components[nc], where, freq_min, freq_max, Nfreq, + chunks = add_dft(components[nc], where, freq, include_dV_and_interp_weights, stored_weight, chunks, sqrt_dV_and_interp_weights,extra_weight,use_centered_grid); - return dft_fields(chunks, freq_min, freq_max, Nfreq, where); + return dft_fields(chunks, freq, where); } /***************************************************************/ @@ -764,7 +835,7 @@ cdouble dft_chunk::process_dft_component(int rank, direction *ds, ivec min_corne ? parent->get_eps(loc) : c_conjugate == Permeability ? parent->get_mu(loc) - : dft[Nomega * (chunk_idx++) + num_freq] / stored_weight); + : dft[omega.size() * (chunk_idx++) + num_freq] / stored_weight); if (include_dV_and_interp_weights) dft_val /= (sqrt_dV_and_interp_weights ? sqrt(w) : w); cdouble mode1val = 0.0, mode2val = 0.0; @@ -1043,7 +1114,7 @@ void fields::output_dft_components(dft_chunk **chunklists, int num_chunklists, v const char *HDF5FileName) { int NumFreqs = 0; for (int nc = 0; nc < num_chunklists && NumFreqs == 0; nc++) - if (chunklists[nc]) NumFreqs = chunklists[nc]->Nomega; + if (chunklists[nc]) NumFreqs = chunklists[nc]->omega.size(); // if the volume has zero thickness in one or more directions, the DFT // grid is two pixels thick in those directions, but we want the HDF5 output diff --git a/src/dft_ldos.cpp b/src/dft_ldos.cpp index add49c9b8..315efd8db 100644 --- a/src/dft_ldos.cpp +++ b/src/dft_ldos.cpp @@ -23,19 +23,20 @@ using namespace std; namespace meep { dft_ldos::dft_ldos(double freq_min, double freq_max, int Nfreq) { - if (Nfreq <= 1) { - omega_min = (freq_min + freq_max) * pi; - domega = 0; - Nomega = 1; - } - else { - omega_min = freq_min * 2 * pi; - domega = (freq_max - freq_min) * 2 * pi / (Nfreq - 1); - Nomega = Nfreq; - } - Fdft = new complex[Nomega]; - Jdft = new complex[Nomega]; - for (int i = 0; i < Nomega; ++i) + freq = meep::linspace(freq_min, freq_max, Nfreq); + Fdft = new complex[Nfreq]; + Jdft = new complex[Nfreq]; + for (int i = 0; i < Nfreq; ++i) + Fdft[i] = Jdft[i] = 0.0; + Jsum = 1.0; +} + +dft_ldos::dft_ldos(const std::vector freq_) { + const int Nfreq = freq_.size(); + freq = freq_; + Fdft = new complex[Nfreq]; + Jdft = new complex[Nfreq]; + for (int i = 0; i < Nfreq; ++i) Fdft[i] = Jdft[i] = 0.0; Jsum = 1.0; } @@ -55,24 +56,27 @@ double *dft_ldos::ldos() const { * -0.5 // power = -1/2 Re[E* J] / (Jsum_all * Jsum_all); // normalize to unit-integral current - double *sum = new double[Nomega]; - for (int i = 0; i < Nomega; ++i) /* 4/pi * work done by unit dipole */ + const int Nfreq = freq.size(); + double *sum = new double[Nfreq]; + for (int i = 0; i < Nfreq; ++i) /* 4/pi * work done by unit dipole */ sum[i] = scale * real(Fdft[i] * conj(Jdft[i])) / abs2(Jdft[i]); - double *out = new double[Nomega]; - sum_to_all(sum, out, Nomega); + double *out = new double[Nfreq]; + sum_to_all(sum, out, Nfreq); delete[] sum; return out; } complex *dft_ldos::F() const { - complex *out = new complex[Nomega]; - sum_to_all(Fdft, out, Nomega); + const int Nfreq = freq.size(); + complex *out = new complex[Nfreq]; + sum_to_all(Fdft, out, Nfreq); return out; } complex *dft_ldos::J() const { - complex *out = new complex[Nomega]; - sum_to_all(Jdft, out, Nomega); + const int Nfreq = freq.size(); + complex *out = new complex[Nfreq]; + sum_to_all(Jdft, out, Nfreq); return out; } @@ -129,9 +133,9 @@ void dft_ldos::update(fields &f) { } } } - for (int i = 0; i < Nomega; ++i) { - complex Ephase = polar(1.0, (omega_min + i * domega) * f.time()) * scale; - complex Hphase = polar(1.0, (omega_min + i * domega) * (f.time() - f.dt / 2)) * scale; + for (int i = 0; i < freq.size(); ++i) { + complex Ephase = polar(1.0, 2*pi * freq[i] * f.time()) * scale; + complex Hphase = polar(1.0, 2*pi * freq[i] * (f.time() - f.dt / 2)) * scale; Fdft[i] += Ephase * EJ + Hphase * HJ; // NOTE: take only 1st time dependence: assumes all sources have same J(t) diff --git a/src/meep.hpp b/src/meep.hpp index 4d95e8289..1f4444b1b 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -62,6 +62,9 @@ class h5file; // Defined in monitor.cpp void matrix_invert(std::complex (&Vinv)[9], std::complex (&V)[9]); +// Defined in dft.cpp +std::vector linspace(double freq_min, double freq_max, int Nfreq); + double pml_quadratic_profile(double, void *); /* generic base class, only used by subclassing: represents susceptibility @@ -371,8 +374,6 @@ class multilevel_susceptibility : public susceptibility { realnum *sigmat; // 5*T transition-specific sigma-diagonal factors }; -class grace; - // h5file.cpp: HDF5 file I/O. Most users, if they use this // class at all, will only use the constructor to open the file, and // will otherwise use the fields::output_hdf5 functions. @@ -1009,8 +1010,7 @@ class dft_chunk { void operator-=(const dft_chunk &chunk); // the frequencies to loop_in_chunks - double omega_min, domega; - int Nomega; + std::vector omega; component c; // component to DFT (possibly transformed by symmetry) @@ -1081,8 +1081,11 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const // dft.cpp (normally created with fields::add_dft_flux) class dft_flux { public: - dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, - double fmax, int Nf, const volume &where_, direction normal_direction_, + dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, + double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, + bool use_symmetry_); + dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, + const std::vector freq_, const volume &where_, direction normal_direction_, bool use_symmetry_); dft_flux(const dft_flux &f); @@ -1103,8 +1106,7 @@ class dft_flux { void remove(); - double freq_min, dfreq; - int Nfreq; + std::vector freq; dft_chunk *E, *H; component cE, cH; volume where; @@ -1115,8 +1117,10 @@ class dft_flux { // dft.cpp (normally created with fields::add_dft_energy) class dft_energy { public: - dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, + dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double freq_min, double freq_max, int Nf, const volume &where_); + dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, const std::vector freq_, + const volume &where_); dft_energy(const dft_energy &f); double *electric(); @@ -1140,8 +1144,7 @@ class dft_energy { void remove(); - double freq_min, dfreq; - int Nfreq; + std::vector freq; dft_chunk *E, *H, *D, *B; volume where; }; @@ -1151,6 +1154,8 @@ class dft_force { public: dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_); + dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, const std::vector freq_, + const volume &where_); dft_force(const dft_force &f); double *force(); @@ -1167,8 +1172,7 @@ class dft_force { void remove(); - double freq_min, dfreq; - int Nfreq; + std::vector freq; dft_chunk *offdiag1, *offdiag2, *diag; volume where; }; @@ -1181,6 +1185,9 @@ class dft_near2far { dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]); + dft_near2far(dft_chunk *F, const std::vector freq_, double eps, double mu, + const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], + const double periodic_k_[2], const double period_[2]); dft_near2far(const dft_near2far &f); /* return an array (Ex,Ey,Ez,Hx,Hy,Hz) x Nfreq of the far fields at x */ @@ -1215,8 +1222,7 @@ class dft_near2far { void remove(); - double freq_min, dfreq; - int Nfreq; + std::vector freq; dft_chunk *F; double eps, mu; volume where; @@ -1234,6 +1240,7 @@ class dft_near2far { class dft_ldos { public: dft_ldos(double freq_min, double freq_max, int Nfreq); + dft_ldos(const std::vector freq); ~dft_ldos() { delete[] Fdft; delete[] Jdft; @@ -1243,27 +1250,25 @@ class dft_ldos { double *ldos() const; // returns array of Nomega values (after last timestep) std::complex *F() const; // returns Fdft std::complex *J() const; // returns Jdft + std::vector freq; private: std::complex *Fdft; // Nomega array of field * J*(x) DFT values std::complex *Jdft; // Nomega array of J(t) DFT values double Jsum; // sum of |J| over all points -public: - double omega_min, domega; - int Nomega; }; // dft.cpp (normally created with fields::add_dft_fields) class dft_fields { public: dft_fields(dft_chunk *chunks, double freq_min, double freq_max, int Nfreq, const volume &where); + dft_fields(dft_chunk *chunks, const std::vector freq_, const volume &where); void scale_dfts(std::complex scale); void remove(); - double freq_min, dfreq; - int Nfreq; + std::vector freq; dft_chunk *chunks; volume where; }; @@ -1727,23 +1732,41 @@ class fields { bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, int vc = 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); dft_chunk *add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq); + dft_chunk *add_dft_pt(component c, const vec &where, const std::vector freq); dft_chunk *add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool include_dV = true); + dft_chunk *add_dft(const volume_list *where, const std::vector freq, + bool include_dV = true); void update_dfts(); dft_flux add_dft_flux(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); + dft_flux add_dft_flux(const volume_list *where, const std::vector freq, + bool use_symmetry = true); dft_flux add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); + dft_flux add_dft_flux(direction d, const volume &where, const std::vector freq, + bool use_symmetry = true); 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); + dft_flux add_dft_flux_plane(const volume &where, const std::vector freq); // 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); + dft_flux add_mode_monitor(direction d, const volume &where, const std::vector freq); 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); + dft_fields add_dft_fields(component *components, int num_components, const volume where, + const std::vector freq, bool use_centered_grid=true); /********************************************************/ /* process_dft_component is an intermediate-level */ @@ -1790,13 +1813,17 @@ class fields { std::complex overlaps[2]); dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq); + dft_energy add_dft_energy(const volume_list *where, const std::vector freq); // stress.cpp dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq); + dft_force add_dft_force(const volume_list *where, const std::vector freq); // near2far.cpp dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, int Nfreq, int Nperiods = 1); + dft_near2far add_dft_near2far(const volume_list *where, const std::vector freq, + int Nperiods = 1); // monitor.cpp std::complex get_chi1inv(component, direction, const vec &loc, double omega = 0, bool parallel = true) const; diff --git a/src/mpb.cpp b/src/mpb.cpp index ea1637801..74ee475ed 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -758,9 +758,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in double *vgrp, kpoint_func user_kpoint_func, void *user_kpoint_data, vec *kpoints, vec *kdom_list, double *cscale, direction d) { - double freq_min = flux.freq_min; - double dfreq = flux.dfreq; - int num_freqs = flux.Nfreq; + int num_freqs = flux.freq.size(); bool match_frequency = true; if (flux.use_symmetry && S.multiplicity() > 1 && parity == 0) @@ -779,12 +777,11 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in /*- call mpb to compute the eigenmode --------------------------*/ /*--------------------------------------------------------------*/ int band_num = bands[nb]; - double freq = freq_min + nf * dfreq; double kdom[3]; - if (user_kpoint_func) kpoint = user_kpoint_func(freq, band_num, user_kpoint_data); + if (user_kpoint_func) kpoint = user_kpoint_func(flux.freq[nf], band_num, user_kpoint_data); am_now_working_on(MPBTime); void *mode_data = - get_eigenmode(freq, d, flux.where, eig_vol, band_num, kpoint, match_frequency, parity, + get_eigenmode(flux.freq[nf], d, flux.where, eig_vol, band_num, kpoint, match_frequency, parity, eig_resolution, eigensolver_tol, kdom, (void **)&mdata); finished_working(); if (!mode_data) { // mode not found, assume evanescent diff --git a/src/near2far.cpp b/src/near2far.cpp index 1628d77d7..347a08267 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -33,10 +33,22 @@ dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, doub const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]) - : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + : F(F_), eps(eps_), mu(mu_), where(where_) { + freq = meep::linspace(fmin, fmax, Nf); + for (int i = 0; i < 2; ++i) { + periodic_d[i] = periodic_d_[i]; + periodic_n[i] = periodic_n_[i]; + periodic_k[i] = periodic_k_[i]; + period[i] = period_[i]; + } +} + +dft_near2far::dft_near2far(dft_chunk *F_, const std::vector freq_, double eps_, double mu_, + const volume &where_, const direction periodic_d_[2], + const int periodic_n_[2], const double periodic_k_[2], + const double period_[2]) + : F(F_), eps(eps_), mu(mu_), where(where_) { + freq = freq_; for (int i = 0; i < 2; ++i) { periodic_d[i] = periodic_d_[i]; periodic_n[i] = periodic_n_[i]; @@ -46,8 +58,8 @@ dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, doub } dft_near2far::dft_near2far(const dft_near2far &f) - : freq_min(f.freq_min), dfreq(f.dfreq), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), - where(f.where) { + : F(f.F), eps(f.eps), mu(f.mu), where(f.where) { + freq = f.freq; for (int i = 0; i < 2; ++i) { periodic_d[i] = f.periodic_d[i]; periodic_n[i] = f.periodic_n[i]; @@ -315,11 +327,12 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { abort("only 2d or 3d or cylindrical far-field computation is supported"); greenfunc green = x.dim == D2 ? green2d : green3d; + const int Nfreq = freq.size(); for (int i = 0; i < 6 * Nfreq; ++i) EH[i] = 0.0; for (dft_chunk *f = F; f; f = f->next_in_dft) { - assert(Nfreq == f->Nomega); + assert(Nfreq == f->omega.size()); component c0 = component(f->vc); /* equivalent source component */ @@ -329,7 +342,6 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { #endif for (int i = 0; i < Nfreq; ++i) { std::complex EH6[6]; - double freq = freq_min + i * dfreq; size_t idx_dft = 0; LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) { IVEC_LOOP_LOC(f->fc->gv, x0); @@ -345,9 +357,9 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { double phase = phase0 + i1 * periodic_k[1]; std::complex cphase = std::polar(1.0, phase); if (x.dim == Dcyl) - greencyl(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i], f->fc->m, 1e-3); + greencyl(EH6, x, freq[i], eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i], f->fc->m, 1e-3); else - green(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]); + green(EH6, x, freq[i], eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]); for (int j = 0; j < 6; ++j) EH[i * 6 + j] += EH6[j] * cphase; } @@ -360,6 +372,7 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { std::complex *dft_near2far::farfield(const vec &x) { std::complex *EH, *EH_local; + const int Nfreq = freq.size(); EH_local = new std::complex[6 * Nfreq]; farfield_lowlevel(EH_local, x); EH = new std::complex[6 * Nfreq]; @@ -388,6 +401,7 @@ realnum *dft_near2far::get_farfields_array(const volume &where, int &rank, size_ } if (where.dim == Dcyl) dirs[2] = P; // otherwise Z is listed twice + const int Nfreq = freq.size(); if (N * Nfreq < 1) return NULL; /* nothing to output */ /* 6 x 2 x N x Nfreq array of fields in row-major order */ @@ -450,6 +464,7 @@ void dft_near2far::save_farfields(const char *fname, const char *prefix, const v realnum *EH = get_farfields_array(where, rank, dims, N, resolution); if (!EH) return; /* nothing to output */ + const int Nfreq = freq.size(); /* frequencies are the last dimension */ if (Nfreq > 1) dims[rank++] = Nfreq; @@ -482,6 +497,7 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution) realnum *EH = get_farfields_array(where, rank, dims, N, resolution); + const int Nfreq = freq.size(); double *F = new double[Nfreq]; std::complex ff_EH[6]; std::complex cE[2], cH[2]; @@ -525,6 +541,11 @@ static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fab dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max, int Nfreq, int Nperiods) { + return add_dft_near2far(where, linspace(freq_min, freq_max, Nfreq), Nperiods); +} + +dft_near2far fields::add_dft_near2far(const volume_list *where, const std::vector freq, + int Nperiods) { dft_chunk *F = 0; /* E and H chunks*/ double eps = 0, mu = 0; volume everywhere = where->v; @@ -598,13 +619,12 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double s = j == 0 ? 1 : -1; /* sign of n x c */ if (is_electric(c)) s = -s; - F = add_dft(c, w->v, freq_min, freq_max, Nfreq, true, s * w->weight, F, false, 1.0, false, - c0); + F = add_dft(c, w->v, freq, true, s * w->weight, F, false, 1.0, false, c0); } } } - return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere, periodic_d, periodic_n, + return dft_near2far(F, freq, eps, mu, everywhere, periodic_d, periodic_n, periodic_k, period); } diff --git a/src/stress.cpp b/src/stress.cpp index 4dc05eb76..a2c1c0eb0 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -27,10 +27,17 @@ namespace meep { dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_) : where(where_) { - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - Nfreq = Nf; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + freq = meep::linspace(fmin, fmax, Nf); + offdiag1 = offdiag1_; + offdiag2 = offdiag2_; + diag = diag_; + // where = new volume(where_.get_min_corner(), where_.get_max_corner()); +} + +dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, + const std::vector freq_, const volume &where_) + : where(where_) { + freq = freq_; offdiag1 = offdiag1_; offdiag2 = offdiag2_; diag = diag_; @@ -38,9 +45,7 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag } dft_force::dft_force(const dft_force &f) : where(f.where) { - freq_min = f.freq_min; - Nfreq = f.Nfreq; - dfreq = f.dfreq; + freq = f.freq; offdiag1 = f.offdiag1; offdiag2 = f.offdiag2; diag = f.diag; @@ -82,6 +87,7 @@ static void stress_sum(int Nfreq, double *F, const dft_chunk *F1, const dft_chun } double *dft_force::force() { + const int Nfreq = freq.size(); double *F = new double[Nfreq]; for (int i = 0; i < Nfreq; ++i) F[i] = 0; @@ -129,11 +135,15 @@ void dft_force::scale_dfts(complex scale) { if (diag) diag->scale_dft(scale); } +dft_force fields::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)); +} + /* 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_, double freq_min, double freq_max, - int Nfreq) { +dft_force fields::add_dft_force(const volume_list *where_, const std::vector freq) { dft_chunk *offdiag1 = 0, *offdiag2 = 0, *diag = 0; volume_list *where = S.reduce(where_); @@ -148,28 +158,28 @@ dft_force fields::add_dft_force(const volume_list *where_, double freq_min, doub if (coordinate_mismatch(gv.dim, fd)) abort("coordinate-type mismatch in add_dft_force"); if (fd != nd) { // off-diagaonal stress-tensor terms - offdiag1 = add_dft(direction_component(Ex, fd), where->v, freq_min, freq_max, Nfreq, true, + offdiag1 = add_dft(direction_component(Ex, fd), where->v, freq, true, where->weight, offdiag1); - offdiag2 = add_dft(direction_component(Ex, nd), where->v, freq_min, freq_max, Nfreq, false, + offdiag2 = add_dft(direction_component(Ex, nd), where->v, freq, false, 1.0, offdiag2); - offdiag1 = add_dft(direction_component(Hx, fd), where->v, freq_min, freq_max, Nfreq, true, + offdiag1 = add_dft(direction_component(Hx, fd), where->v, freq, true, where->weight, offdiag1); - offdiag2 = add_dft(direction_component(Hx, nd), where->v, freq_min, freq_max, Nfreq, false, + offdiag2 = add_dft(direction_component(Hx, nd), where->v, freq, false, 1.0, offdiag2); } 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_min, freq_max, Nfreq, true, 1.0, + diag = add_dft(direction_component(Ex, d), where->v, freq, true, 1.0, diag, true, weight1, false); - diag = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, + diag = add_dft(direction_component(Hx, d), where->v, freq, true, 1.0, diag, true, weight1, false); } everywhere = everywhere | where->v; } delete where_save; - return dft_force(offdiag1, offdiag2, diag, freq_min, freq_max, Nfreq, everywhere); + return dft_force(offdiag1, offdiag2, diag, freq, everywhere); } } // namespace meep diff --git a/tests/flux.cpp b/tests/flux.cpp index f52836546..f5cbc6fdd 100644 --- a/tests/flux.cpp +++ b/tests/flux.cpp @@ -180,10 +180,12 @@ int flux_2d(const double xmax, const double ymax, double eps(const vec &)) { around the source...should be positive and equal */ volume box1(vec(xmax / 6 - 0.4, ymax / 6 - 0.2), vec(xmax / 6 + 0.6, ymax / 6 + 0.8)); volume box2(vec(xmax / 6 - 0.9, ymax / 6 - 0.7), vec(xmax / 6 + 1.1, ymax / 6 + 1.3)); - double fmin = 0.23, fmax = 0.27; int Nfreq = 10; - dft_flux flux1 = f.add_dft_flux_box(box1, fmin, fmax, Nfreq); - dft_flux flux2 = f.add_dft_flux_box(box2, fmin, fmax, Nfreq); + double freq_array[] = {0.230, 0.232, 0.238, 0.241, 0.248, 0.254, 0.256, 0.265, 0.269, 0.270}; + // workaround for C++98 which does not support list initialization + const std::vector freq(freq_array, freq_array + sizeof(freq_array)/sizeof(double)); + dft_flux flux1 = f.add_dft_flux_box(box1, freq); + dft_flux flux2 = f.add_dft_flux_box(box2, freq); const double ttot = 130; @@ -213,7 +215,7 @@ int flux_2d(const double xmax, const double ymax, double eps(const vec &)) { double *fl1 = flux1.flux(); double *fl2 = flux2.flux(); for (int i = 0; i < Nfreq; ++i) { - master_printf(" flux(%g) = %g vs. %g (rat. = %g)\n", fmin + i * flux1.dfreq, fl1[i], fl2[i], + master_printf(" flux(%g) = %g vs. %g (rat. = %g)\n", flux1.freq[i], fl1[i], fl2[i], fl1[i] / fl2[i]); if (!compare(fl1[i], fl2[i], 0.09, 0, "Flux spectrum")) return 0; } @@ -281,7 +283,7 @@ int flux_cyl(const double rmax, const double zmax, double eps(const vec &), int double *fl1 = flux1.flux(); double *fl2 = flux2.flux(); for (int i = 0; i < Nfreq; ++i) { - master_printf(" flux(%g) = %g vs. %g (rat. = %g)\n", fmin + i * flux1.dfreq, fl1[i], fl2[i], + master_printf(" flux(%g) = %g vs. %g (rat. = %g)\n", flux1.freq[i], fl1[i], fl2[i], fl1[i] / fl2[i]); if (!compare(fl1[i], fl2[i], 0.08, 0, "Flux spectrum")) return 0; } diff --git a/tests/near2far.cpp b/tests/near2far.cpp index 4c105bad9..b1154d377 100644 --- a/tests/near2far.cpp +++ b/tests/near2far.cpp @@ -59,7 +59,6 @@ int check_cyl(double sr, double sz, double a) { greencyl(EH, x, w, 2.0, 1.0, x0, c0, 1.0, m, 1e-6); F0[i] = EH[EHcomp[c]] * 2.0*pi*x0.r(); // Ey = Ep for \phi = 0 (rz = xz) plane double d = abs(F0[i] - F[i]); - double f = abs(F[i]); double f0 = abs(F0[i]); diff += d * d; dot0 += f0 * f0;