From 15d17674611e73a40649f69befd831c73aaf3e56 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Mon, 6 Apr 2020 07:23:33 -0700 Subject: [PATCH] rename omega parameter to frequency (#1171) * rename omega parameter to frequency * fix --- doc/docs/FAQ.md | 2 +- doc/docs/Python_User_Interface.md | 28 ++++++------- python/meep.i | 16 ++++---- python/simulation.py | 52 +++++++++++++++++------- python/tests/dispersive_eigenmode.py | 18 ++++----- python/tests/mode_coeffs.py | 2 +- python/tests/mpb.py | 6 +-- python/visualization.py | 18 ++++----- src/array_slice.cpp | 44 ++++++++++---------- src/h5fields.cpp | 52 ++++++++++++------------ src/meep.hpp | 60 ++++++++++++++-------------- src/monitor.cpp | 48 +++++++++++----------- src/mpb.cpp | 52 ++++++++++++------------ 13 files changed, 210 insertions(+), 188 deletions(-) diff --git a/doc/docs/FAQ.md b/doc/docs/FAQ.md index 14104efb5..837586612 100644 --- a/doc/docs/FAQ.md +++ b/doc/docs/FAQ.md @@ -407,7 +407,7 @@ Note: a simple approach to reduce the cost of the DTFT computation is to reduce ### Does Meep support shared-memory parallelism? -You can always run the MPI parallel Meep on a shared-memory machine, and some MPI implementations take special advantage of shared memory communications. Meep currently also provides limited support for [multithreading](https://en.wikipedia.org/wiki/Thread_(computing)#Multithreading) via OpenMP on a single, shared-memory, multi-core machine to speed up *multi-frequency* [near-to-far field](Python_User_Interface.md#near-to-far-field-spectra) calculations involving `get_farfields` or `output_farfields`. +Yes. Meep provides limited support for [multithreading](https://en.wikipedia.org/wiki/Thread_(computing)#Multithreading) via OpenMP on a single, shared-memory, multi-core machine to speed up *multi-frequency* [near-to-far field](Python_User_Interface.md#near-to-far-field-spectra) calculations involving `get_farfields` or `output_farfields`. Also, you can always run the MPI parallel Meep on a shared-memory machine, and some MPI implementations take special advantage of shared memory communications. ### Why does the time-stepping rate fluctuate erratically for jobs running on a shared-memory system? diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 6d1ee296b..4ce844dd4 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -1042,9 +1042,9 @@ Sets the condition of the boundary on the specified side in the specified direct — Given a `component` or `derived_component` constant `c` and a `Vector3` `pt`, returns the value of that component at that point. -**`get_epsilon_point(pt, omega=0)`** +**`get_epsilon_point(pt, frequency=0)`** — -Given a frequency `omega` and a `Vector3` `pt`, returns the average eigenvalue of the permittivity tensor at that location and frequency. If `omega` is non-zero, the result is complex valued; otherwise it is the real, frequency-independent part of ε (the $\omega\to\infty$ limit). +Given a frequency `frequency` and a `Vector3` `pt`, returns the average eigenvalue of the permittivity tensor at that location and frequency. If `frequency` is non-zero, the result is complex valued; otherwise it is the real, frequency-independent part of ε (the $\omega\to\infty$ limit). **`initialize_field(c, func)`** — @@ -1220,7 +1220,7 @@ Given a flux object and list of band indices, return a `namedtuple` with the fol + `kdom`: a list of `mp.Vector3`s of the mode's dominant wavevector. + `cscale`: a NumPy array of each mode's scaling coefficient. Useful for adjoint calculations. -The flux object should be created using `add_mode_monitor`. (You could also use `add_flux`, but with `add_flux` you need to be more careful about symmetries that bisect the flux plane: the `add_flux` object should only be used with `get_eigenmode_coefficients` for modes of the same symmetry, e.g. constrained via `eig_parity`. On the other hand, the performance of `add_flux` planes benefits more from symmetry.) `eig_vol` is the volume passed to [MPB](https://mpb.readthedocs.io) for the eigenmode calculation (based on interpolating the discretized materials from the Yee grid); in most cases this will simply be the volume over which the frequency-domain fields are tabulated, which is the default (i.e. `flux.where`). `eig_parity` should be one of [`mp.NO_PARITY` (default), `mp.EVEN_Z`, `mp.ODD_Z`, `mp.EVEN_Y`, `mp.ODD_Y`]. It is the parity (= polarization in 2d) of the mode to calculate, assuming the structure has $z$ and/or $y$ mirror symmetry *in the source region*, just as for `EigenmodeSource` above. If the structure has both $y$ and $z$ mirror symmetry, you can combine more than one of these, e.g. `EVEN_Z+ODD_Y`. Default is `NO_PARITY`, in which case MPB computes all of the bands which will still be even or odd if the structure has mirror symmetry, of course. This is especially useful in 2d simulations to restrict yourself to a desired polarization. `eig_resolution` is the spatial resolution to use in MPB for the eigenmode calculations. This defaults to twice the Meep `resolution` in which case the structure is linearly interpolated from the Meep pixels. `eig_tolerance` is the tolerance to use in the MPB eigensolver. MPB terminates when the eigenvalues stop changing to less than this fractional tolerance. Defaults to `1e-12`. (Note that this is the tolerance for the frequency eigenvalue ω; the tolerance for the mode profile is effectively the square root of this.) For examples, see [Tutorial/Mode Decomposition](Python_Tutorials/Mode_Decomposition.md). +The flux object should be created using `add_mode_monitor`. (You could also use `add_flux`, but with `add_flux` you need to be more careful about symmetries that bisect the flux plane: the `add_flux` object should only be used with `get_eigenmode_coefficients` for modes of the same symmetry, e.g. constrained via `eig_parity`. On the other hand, the performance of `add_flux` planes benefits more from symmetry.) `eig_vol` is the volume passed to [MPB](https://mpb.readthedocs.io) for the eigenmode calculation (based on interpolating the discretized materials from the Yee grid); in most cases this will simply be the volume over which the frequency-domain fields are tabulated, which is the default (i.e. `flux.where`). `eig_parity` should be one of [`mp.NO_PARITY` (default), `mp.EVEN_Z`, `mp.ODD_Z`, `mp.EVEN_Y`, `mp.ODD_Y`]. It is the parity (= polarization in 2d) of the mode to calculate, assuming the structure has $z$ and/or $y$ mirror symmetry *in the source region*, just as for `EigenModeSource` above. If the structure has both $y$ and $z$ mirror symmetry, you can combine more than one of these, e.g. `EVEN_Z+ODD_Y`. Default is `NO_PARITY`, in which case MPB computes all of the bands which will still be even or odd if the structure has mirror symmetry, of course. This is especially useful in 2d simulations to restrict yourself to a desired polarization. `eig_resolution` is the spatial resolution to use in MPB for the eigenmode calculations. This defaults to twice the Meep `resolution` in which case the structure is linearly interpolated from the Meep pixels. `eig_tolerance` is the tolerance to use in the MPB eigensolver. MPB terminates when the eigenvalues stop changing to less than this fractional tolerance. Defaults to `1e-12`. (Note that this is the tolerance for the frequency eigenvalue ω; the tolerance for the mode profile is effectively the square root of this.) For examples, see [Tutorial/Mode Decomposition](Python_Tutorials/Mode_Decomposition.md). Technically, MPB computes `ωₙ(k)` and then inverts it with Newton's method to find the wavevector `k` normal to `eig_vol` and mode for a given frequency; in rare cases (primarily waveguides with *nonmonotonic* dispersion relations, which doesn't usually happen in simple dielectric waveguides), MPB may need you to supply an initial "guess" for `k` in order for this Newton iteration to converge. You can supply this initial guess with `kpoint_func`, which is a function `kpoint_func(f, n)` that supplies a rough initial guess for the `k` of band number `n` at frequency `f = ω/2π`. (By default, the **k** components in the plane of the `eig_vol` region are zero. However, if this region spans the *entire* cell in some directions, and the cell has Bloch-periodic boundary conditions via the `k_point` parameter, then the mode's **k** components in those directions will match `k_point` so that the mode satisfies the Meep boundary conditions, regardless of `kpoint_func`.) If `direction` is set to `mp.NO_DIRECTION`, then `kpoint_func` is not only the initial guess and the search direction of the **k** vectors, but is also taken to be the direction of the waveguide, allowing you to [detect modes in oblique waveguides](Python_Tutorials/Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide) (not perpendicular to the flux plane). @@ -1232,12 +1232,12 @@ Similar to `add_flux`, but for use with `get_eigenmode_coefficients`. `add_mode_monitor` works properly with arbitrary symmetries, but may be suboptimal because the Fourier-transformed region does not exploit the symmetry. As an optimization, if you have a mirror plane that bisects the mode monitor, you can instead use `add_flux` to gain a factor of two, but in that case you *must* also pass the corresponding `eig_parity` to `get_eigenmode_coefficients` in order to only compute eigenmodes with the corresponding mirror symmetry. -**`get_eigenmode(freq, direction, where, band_num, kpoint, eig_vol=None, match_frequency=True, parity=mp.NO_PARITY, resolution=0, eigensolver_tol=1e-12)`** +**`get_eigenmode(frequency, direction, where, band_num, kpoint, eig_vol=None, match_frequency=True, parity=mp.NO_PARITY, resolution=0, eigensolver_tol=1e-12)`** — The parameters of this routine are the same as that of `get_eigenmode_coefficients` or `EigenModeSource`, but this function returns an object that can be used to inspect the computed mode. In particular, it returns an `EigenmodeData` instance with the following fields: + `band_num`: same as the `band_num` parameter -+ `freq`: the computed frequency, same as the `freq` input parameter if `match_frequency=True` ++ `frequency`: the computed frequency, same as the `frequency` input parameter if `match_frequency=True` + `group_velocity`: the group velocity of the mode in `direction` + `k`: the Bloch wavevector of the mode in `direction` + `kdom`: the dominant planewave of mode `band_num` @@ -1615,7 +1615,7 @@ fr = mp.FluxRegion(volume=mp.GDSII_vol(fname, layer, zmin, zmax)) This module provides basic visualization functionality for the simulation domain. The spirit of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. -**`Simulation.plot2D(ax=None, output_plane=None, fields=None, labels=False, eps_parameters=None, boundary_parameters=None, source_parameters=None, monitor_parameters=None, field_parameters=None, omega=None)`** +**`Simulation.plot2D(ax=None, output_plane=None, fields=None, labels=False, eps_parameters=None, boundary_parameters=None, source_parameters=None, monitor_parameters=None, field_parameters=None, frequency=None)`** — Plots a 2D cross section of the simulation domain using `matplotlib`. The plot includes the geometry, boundary layers, sources, and monitors. Fields can also be superimposed on a 2D slice. Requires [matplotlib](https://matplotlib.org). Calling this function would look something like: @@ -1670,7 +1670,7 @@ plt.savefig('sim_domain.png') - `cmap='RdBu'`: color map for field pixels - `alpha=0.6`: transparency of fields - `post_process=np.real`: post processing function to apply to fields (must be a function object) -* `omega`: for materials with a [frequency-dependent permittivity](Materials.md#material-dispersion) $\varepsilon(\omega)$, specifies the frequency $\omega$ (in Meep units) of the real part of the permittivity to use in the plot. Defaults to the `frequency` parameter of the [Source](#source) object. +* `frequency`: for materials with a [frequency-dependent permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the frequency $f$ (in Meep units) of the real part of the permittivity to use in the plot. Defaults to the `frequency` parameter of the [Source](#source) object. **`Simulation.plot3D()`** — Uses Mayavi to render a 3D simulation domain. The simulation object must be 3D. Can also be embedded in Jupyter notebooks. @@ -1797,13 +1797,13 @@ The most common step function is an output function, which outputs some field co Note that although the various field components are stored at different places in the [Yee lattice](Yee_Lattice.md), when they are outputted they are all linearly interpolated to the same grid: to the points at the *centers* of the Yee cells, i.e. $(i+0.5,j+0.5,k+0.5)\cdotΔ$ in 3d. -**`output_epsilon(omega=0)`** +**`output_epsilon(frequency=0)`** — -Given a frequency `omega`, output ε (relative permittivity); for an anisotropic ε tensor the output is the [harmonic mean](https://en.wikipedia.org/wiki/Harmonic_mean) of the ε eigenvalues. If `omega` is non-zero, the output is complex; otherwise it is the real, frequency-independent part of ε (the $\omega\to\infty$ limit). +Given a frequency `frequency`, output ε (relative permittivity); for an anisotropic ε tensor the output is the [harmonic mean](https://en.wikipedia.org/wiki/Harmonic_mean) of the ε eigenvalues. If `frequency` is non-zero, the output is complex; otherwise it is the real, frequency-independent part of ε (the $\omega\to\infty$ limit). -**`output_mu(omega=0)`** +**`output_mu(frequency=0)`** — -Given a frequency `omega`, output μ (relative permeability); for an anisotropic μ tensor the output is the [harmonic mean](https://en.wikipedia.org/wiki/Harmonic_mean) of the μ eigenvalues. If `omega` is non-zero, the output is complex; otherwise it is the real, frequency-independent part of μ (the $\omega\to\infty$ limit). +Given a frequency `frequency`, output μ (relative permeability); for an anisotropic μ tensor the output is the [harmonic mean](https://en.wikipedia.org/wiki/Harmonic_mean) of the μ eigenvalues. If `frequency` is non-zero, the output is complex; otherwise it is the real, frequency-independent part of μ (the $\omega\to\infty$ limit). **`Simulation.output_dft(dft_fields, fname)`** — @@ -1854,7 +1854,7 @@ See also [Field Functions](Field_Functions.md), and [Synchronizing the Magnetic The output functions described above write the data for the fields and materials for the entire cell to an HDF5 file. This is useful for post-processing as you can later read in the HDF5 file to obtain field/material data as a NumPy array. However, in some cases it is convenient to bypass the disk altogether to obtain the data *directly* in the form of a NumPy array without writing/reading HDF5 files. Additionally, you may want the field/material data on just a subregion (or slice) of the entire volume. This functionality is provided by the `get_array` method which takes as input a subregion of the cell and the field/material component. The method returns a NumPy array containing values of the field/material at the current simulation time. ```python - get_array(vol=None, center=None, size=None, component=mp.Ez, cmplx=False, arr=None, omega=0) + get_array(vol=None, center=None, size=None, component=mp.Ez, cmplx=False, arr=None, frequency=0) ``` with the following input parameters: @@ -1869,9 +1869,9 @@ with the following input parameters: + `arr`: optional field to pass a pre-allocated NumPy array of the correct size, which will be overwritten with the field/material data instead of allocating a new array. Normally, this will be the array returned from a previous call to `get_array` for a similar slice, allowing one to re-use `arr` (e.g., when fetching the same slice repeatedly at different times). -+ `omega`: optional frequency point over which the average eigenvalue of the dielectric and permeability tensors are evaluated (defaults to 0). ++ `frequency`: optional frequency point over which the average eigenvalue of the dielectric and permeability tensors are evaluated (defaults to 0). -For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` accept the optional `omega` parameter (defaults to 0). +For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` accept the optional `frequency` parameter (defaults to 0). **Note on array-slice dimensions:** The routines `get_epsilon`, `get_Xfield_z`, etc. use as default `size=meep.Simulation.fields.total_volume()` which for simulations involving Bloch-periodic boundaries (via `k_point`) will result in arrays that have slightly *different* dimensions than e.g. `get_array(center=meep.Vector3(), size=cell_size, component=meep.Dielectric`, etc. (i.e., the slice spans the entire cell volume `cell_size`). Neither of these approaches is "wrong", they are just slightly different methods of fetching the boundaries. The key point is that if you pass the same value for the `size` parameter, or use the default, the slicing routines always give you the same-size array for all components. You should *not* try to predict the exact size of these arrays; rather, you should simply rely on Meep's output. diff --git a/python/meep.i b/python/meep.i index a99ccc631..f3c5d0f49 100644 --- a/python/meep.i +++ b/python/meep.i @@ -53,7 +53,7 @@ namespace meep { vec center; amplitude_function amp_func; int band_num; - double omega; + double frequency; double group_velocity; }; } @@ -519,12 +519,12 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where } #ifdef HAVE_MPB -meep::eigenmode_data *_get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const meep::volume where, +meep::eigenmode_data *_get_eigenmode(meep::fields *f, double frequency, meep::direction d, const meep::volume where, const meep::volume eig_vol, int band_num, const meep::vec &_kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, double kdom[3]) { - void *data = f->get_eigenmode(omega_src, d, where, eig_vol, band_num, _kpoint, match_frequency, + void *data = f->get_eigenmode(frequency, d, where, eig_vol, band_num, _kpoint, match_frequency, parity, resolution, eigensolver_tol, kdom); return (meep::eigenmode_data *)data; } @@ -538,11 +538,11 @@ PyObject *_get_eigenmode_Gk(meep::eigenmode_data *emdata) { } #else -void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const meep::volume where, +void _get_eigenmode(meep::fields *f, double frequency, meep::direction d, const meep::volume where, const meep::volume eig_vol, int band_num, const meep::vec &_kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, double kdom[3]) { - (void) f; (void) omega_src; (void) d; (void) where; (void) eig_vol; (void) band_num; (void) _kpoint; + (void) f; (void) frequency; (void) d; (void) where; (void) eig_vol; (void) band_num; (void) _kpoint; (void) match_frequency; (void) parity; (void) resolution; (void) eigensolver_tol; (void) kdom; meep::abort("Must compile Meep with MPB for get_eigenmode"); @@ -1328,12 +1328,12 @@ struct eigenmode_data { vec center; amplitude_function amp_func; int band_num; - double omega; + double frequency; double group_velocity; }; } -meep::eigenmode_data *_get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const meep::volume where, +meep::eigenmode_data *_get_eigenmode(meep::fields *f, double frequency, meep::direction d, const meep::volume where, const meep::volume eig_vol, int band_num, const meep::vec &_kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, double kdom[3]); @@ -1346,7 +1346,7 @@ PyObject *_get_eigenmode_Gk(meep::eigenmode_data *emdata); } #else -void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const meep::volume where, +void _get_eigenmode(meep::fields *f, double frequency, meep::direction d, const meep::volume where, const meep::volume eig_vol, int band_num, const meep::vec &_kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, double kdom[3]); diff --git a/python/simulation.py b/python/simulation.py index 3aa3b9c75..6aedfb873 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1427,9 +1427,12 @@ def get_field_point(self, c, pt): v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical) return self.fields.get_field_from_comp(c, v3) - def get_epsilon_point(self, pt, omega = 0): + def get_epsilon_point(self, pt, frequency=0, omega=0): v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical) - return self.fields.get_eps(v3,omega) + if omega != 0: + frequency = omega + warnings.warn("get_epsilon_point: omega has been deprecated; use frequency instead", RuntimeWarning) + return self.fields.get_eps(v3,frequency) def get_filename_prefix(self): if isinstance(self.filename_prefix, str): @@ -2006,7 +2009,7 @@ def _add_fluxish_stuff(self, add_dft_stuff, freq, stufflist, *args): return stuff - def output_component(self, c, h5file=None, omega=0): + def output_component(self, c, h5file=None, frequency=0, omega=0): if self.fields is None: raise RuntimeError("Fields must be initialized before calling output_component") @@ -2014,7 +2017,11 @@ def output_component(self, c, h5file=None, omega=0): h5 = self.output_append_h5 if h5file is None else h5file append = h5file is None and self.output_append_h5 is not None - self.fields.output_hdf5(c, vol, h5, append, self.output_single_precision,self.get_filename_prefix(), omega) + if omega != 0: + frequency = omega + warnings.warn("output_component: omega has been deprecated; use frequency instead", RuntimeWarning) + + self.fields.output_hdf5(c, vol, h5, append, self.output_single_precision,self.get_filename_prefix(), frequency) if h5file is None: nm = self.fields.h5file_name(mp.component_name(c), self.get_filename_prefix(), True) @@ -2044,7 +2051,7 @@ def h5topng(self, rm_h5, option, *step_funcs): cmd = re.sub(r'\$EPS', self.last_eps_filename, opts) return convert_h5(rm_h5, cmd, *step_funcs) - def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None, arr=None, omega = 0): + def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None, arr=None, frequency=0, omega=0): if component is None: raise ValueError("component is required") if isinstance(component, mp.Volume) or isinstance(component, mp.volume): @@ -2061,8 +2068,12 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None dims = [s for s in dim_sizes if s != 0] + if omega != 0: + frequency = omega + warnings.warn("get_array: omega has been deprecated; use frequency instead", RuntimeWarning) + if cmplx is None: - cmplx = omega != 0 or (component < mp.Dielectric and not self.fields.is_real) + cmplx = frequency != 0 or (component < mp.Dielectric and not self.fields.is_real) if arr is not None: if cmplx and not np.iscomplexobj(arr): @@ -2079,9 +2090,9 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None arr = np.zeros(dims, dtype=np.complex128 if cmplx else np.float64) if np.iscomplexobj(arr): - self.fields.get_complex_array_slice(v, component, arr, omega) + self.fields.get_complex_array_slice(v, component, arr, frequency) else: - self.fields.get_array_slice(v, component, arr, omega) + self.fields.get_array_slice(v, component, arr, frequency) return arr @@ -2177,7 +2188,7 @@ def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, eig_v 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, + def get_eigenmode(self, frequency, direction, where, band_num, kpoint, eig_vol=None, match_frequency=True, parity=mp.NO_PARITY, resolution=0, eigensolver_tol=1e-12): if self.fields is None: @@ -2191,11 +2202,11 @@ def get_eigenmode(self, freq, direction, where, band_num, kpoint, eig_vol=None, swig_kpoint = mp.vec(kpoint.x, kpoint.y, kpoint.z) kdom = np.zeros(3) - emdata = mp._get_eigenmode(self.fields, freq, direction, where, eig_vol, band_num, swig_kpoint, + emdata = mp._get_eigenmode(self.fields, frequency, direction, where, eig_vol, band_num, swig_kpoint, match_frequency, parity, resolution, eigensolver_tol, kdom) Gk = mp._get_eigenmode_Gk(emdata) - return EigenmodeData(emdata.band_num, emdata.omega, emdata.group_velocity, Gk, + return EigenmodeData(emdata.band_num, emdata.frequency, emdata.group_velocity, Gk, emdata, mp.Vector3(kdom[0], kdom[1], kdom[2])) def output_field_function(self, name, cs, func, real_only=False, h5file=None): @@ -2290,8 +2301,11 @@ def print_times(self): if self.fields: self.fields.print_times() - def get_epsilon(self,omega=0): - return self.get_array(component=mp.Dielectric,omega=omega) + def get_epsilon(self,frequency=0,omega=0): + if omega != 0: + frequency = omega + warnings.warn("get_epsilon: omega has been deprecated; use frequency instead", RuntimeWarning) + return self.get_array(component=mp.Dielectric,frequency=frequency) def get_mu(self): return self.get_array(component=mp.Permeability) @@ -2827,13 +2841,21 @@ def _output_png(sim, todo): def output_epsilon(sim,*step_func_args,**kwargs): + frequency = kwargs.pop('frequency', 0.0) omega = kwargs.pop('omega', 0.0) - sim.output_component(mp.Dielectric,omega=omega) + if omega != 0: + frequency = omega + warnings.warn("output_epsilon: omega has been deprecated; use frequency instead", RuntimeWarning) + sim.output_component(mp.Dielectric,frequency=frequency) def output_mu(sim,*step_func_args,**kwargs): + frequency = kwargs.pop('frequency', 0.0) omega = kwargs.pop('omega', 0.0) - sim.output_component(mp.Permeability,omega=omega) + if omega != 0: + frequency = omega + warnings.warn("output_mu: omega has been deprecated; use frequency instead", RuntimeWarning) + sim.output_component(mp.Permeability,frequency=frequency) def output_hpwr(sim): diff --git a/python/tests/dispersive_eigenmode.py b/python/tests/dispersive_eigenmode.py index af755e0c7..60521b756 100644 --- a/python/tests/dispersive_eigenmode.py +++ b/python/tests/dispersive_eigenmode.py @@ -20,7 +20,7 @@ class TestDispersiveEigenmode(unittest.TestCase): # ----------- Helper Functions ------------ # # ----------------------------------------- # # Directly cals the C++ chi1 routine - def call_chi1(self,material,omega): + def call_chi1(self,material,frequency): sim = mp.Simulation(cell_size=mp.Vector3(1,1,1), default_material=material, @@ -31,10 +31,10 @@ def call_chi1(self,material,omega): chi1inv = np.zeros((3,3),dtype=np.complex128) for i, com in enumerate([mp.Ex,mp.Ey,mp.Ez]): for k, dir in enumerate([mp.X,mp.Y,mp.Z]): - chi1inv[i,k] = sim.structure.get_chi1inv(com,dir,v3,omega) + chi1inv[i,k] = sim.structure.get_chi1inv(com,dir,v3,frequency) n = np.real(np.sqrt(np.linalg.inv(chi1inv.astype(np.complex128)))) - n_actual = np.real(np.sqrt(material.epsilon(omega).astype(np.complex128))) + n_actual = np.real(np.sqrt(material.epsilon(frequency).astype(np.complex128))) np.testing.assert_allclose(n,n_actual) @@ -46,9 +46,9 @@ def setUpClass(cls): def tearDownClass(cls): mp.delete_directory(cls.temp_dir) - def verify_output_and_slice(self,material,omega): + def verify_output_and_slice(self,material,frequency): # Since the slice routines average the diagonals, we need to do that too: - chi1 = material.epsilon(omega).astype(np.complex128) + chi1 = material.epsilon(frequency).astype(np.complex128) chi1inv = np.linalg.inv(chi1) chi1inv = np.diag(chi1inv) N = chi1inv.size @@ -62,13 +62,13 @@ def verify_output_and_slice(self,material,omega): sim.use_output_directory(self.temp_dir) sim.init_sim() - # Check to make sure the get_slice routine is working with omega - n_slice = np.sqrt(np.max(sim.get_epsilon(omega))) + # Check to make sure the get_slice routine is working with frequency + n_slice = np.sqrt(np.max(sim.get_epsilon(frequency))) self.assertAlmostEqual(n,n_slice, places=4) - # Check to make sure h5 output is working with omega + # Check to make sure h5 output is working with frequency filename = os.path.join(self.temp_dir, 'dispersive_eigenmode-eps-000000.00.h5') - mp.output_epsilon(sim,omega=omega) + mp.output_epsilon(sim,frequency=frequency) n_h5 = 0 mp.all_wait() with h5py.File(filename, 'r') as f: diff --git a/python/tests/mode_coeffs.py b/python/tests/mode_coeffs.py index d6b3f61db..d3b650b1e 100644 --- a/python/tests/mode_coeffs.py +++ b/python/tests/mode_coeffs.py @@ -55,7 +55,7 @@ def run_mode_coeffs(self, mode_num, kpoint_func, nf=1, resolution=15): mode_flux = sim.add_flux(fcen, df, nf, mp.FluxRegion(center=mp.Vector3(xm,0), size=mp.Vector3(0,sy-2*dpml))) # sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(-0.5*sx+dpml,0), 1e-10)) - sim.run(until_after_sources=200) + sim.run(until_after_sources=100) ################################################## # If the number of analysis frequencies is >1, we diff --git a/python/tests/mpb.py b/python/tests/mpb.py index c3eed513a..3a315cc01 100644 --- a/python/tests/mpb.py +++ b/python/tests/mpb.py @@ -890,10 +890,10 @@ def test_strip(self): for e, r in zip(expected_z_parities, z_parities): self.assertAlmostEqual(e, r, places=3) - omega = 1 / 1.55 + frequency = 1 / 1.55 - kvals = ms.find_k(mp.NO_PARITY, omega, 1, ms.num_bands, mp.Vector3(1), 1e-3, - omega * 3.45, omega * 0.1, omega * 4, mpb.output_poynting_x, + kvals = ms.find_k(mp.NO_PARITY, frequency, 1, ms.num_bands, mp.Vector3(1), 1e-3, + frequency * 3.45, frequency * 0.1, frequency * 4, mpb.output_poynting_x, mpb.display_yparities, mpb.display_group_velocities) expected_kvals = [ diff --git a/python/visualization.py b/python/visualization.py index 378d2732c..00cbd840c 100644 --- a/python/visualization.py +++ b/python/visualization.py @@ -303,7 +303,7 @@ def sort_points(xy): return ax return ax -def plot_eps(sim,ax,output_plane=None,eps_parameters=None,omega=0): +def plot_eps(sim,ax,output_plane=None,eps_parameters=None,frequency=0): if sim.structure is None: sim.init_sim() @@ -342,7 +342,7 @@ def plot_eps(sim,ax,output_plane=None,eps_parameters=None,omega=0): else: raise ValueError("A 2D plane has not been specified...") - eps_data = np.rot90(np.real(sim.get_array(center=center, size=cell_size, component=mp.Dielectric, omega=omega))) + eps_data = np.rot90(np.real(sim.get_array(center=center, size=cell_size, component=mp.Dielectric, frequency=frequency))) if mp.am_master(): ax.imshow(eps_data, extent=extent, **eps_parameters) ax.set_xlabel(xlabel) @@ -514,7 +514,7 @@ def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None) def plot2D(sim,ax=None, output_plane=None, fields=None, labels=False, eps_parameters=None,boundary_parameters=None, source_parameters=None,monitor_parameters=None, - field_parameters=None, omega=None): + field_parameters=None, frequency=None): # Initialize the simulation if sim.structure is None: @@ -524,21 +524,21 @@ def plot2D(sim,ax=None, output_plane=None, fields=None, labels=False, from matplotlib import pyplot as plt ax = plt.gca() # Determine a frequency to plot all epsilon - if omega is None: + if frequency is None: try: # Continuous sources - omega = sim.sources[0].frequency + frequency = sim.sources[0].frequency except: try: # Gaussian sources - omega = sim.sources[0].src.frequency + frequency = sim.sources[0].src.frequency except: try: # Custom sources - omega = sim.sources[0].src.center_frequency + frequency = sim.sources[0].src.center_frequency except: # No sources - omega = 0 + frequency = 0 # validate the output plane to ensure proper 2D coordinates from meep.simulation import Volume @@ -546,7 +546,7 @@ def plot2D(sim,ax=None, output_plane=None, fields=None, labels=False, output_plane = Volume(center=sim_center,size=sim_size) # Plot geometry - ax = plot_eps(sim,ax,output_plane=output_plane,eps_parameters=eps_parameters,omega=omega) + ax = plot_eps(sim,ax,output_plane=output_plane,eps_parameters=eps_parameters,frequency=frequency) # Plot boundaries ax = plot_boundaries(sim,ax,output_plane=output_plane,boundary_parameters=boundary_parameters) diff --git a/src/array_slice.cpp b/src/array_slice.cpp index 1aed17937..f5920d596 100644 --- a/src/array_slice.cpp +++ b/src/array_slice.cpp @@ -66,7 +66,7 @@ typedef struct { cdouble *fields; ptrdiff_t *offsets; - double omega; + double frequency; int ninveps; component inveps_cs[3]; @@ -282,7 +282,7 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr ptrdiff_t *off = data->offsets; component *cS = data->cS; - double omega = data->omega; + double frequency = data->frequency; complex *fields = data->fields, *ph = data->ph; const component *iecs = data->inveps_cs; const direction *ieds = data->inveps_ds; @@ -325,11 +325,11 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr else if (cS[i] == Dielectric) { complex tr(0.0,0.0); for (int k = 0; k < data->ninveps; ++k) { - tr += (fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx, omega) + - fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k], omega) + - fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[1 + 2 * k], omega) + + tr += (fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx, frequency) + + fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k], frequency) + + fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[1 + 2 * k], frequency) + fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k] + ieos[1 + 2 * k], - omega)); + frequency)); if (abs(tr) == 0.0) tr += 4.0; // default inveps == 1 } fields[i] = (4.0 * data->ninveps) / tr; @@ -337,11 +337,11 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr else if (cS[i] == Permeability) { complex tr(0.0,0.0); for (int k = 0; k < data->ninvmu; ++k) { - tr += (fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx, omega) + - fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k], omega) + - fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[1 + 2 * k], omega) + + tr += (fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx, frequency) + + fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k], frequency) + + fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[1 + 2 * k], frequency) + fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k] + imos[1 + 2 * k], - omega)); + frequency)); if (abs(tr) == 0.0) tr += 4.0; // default invmu == 1 } fields[i] = (4.0 * data->ninvmu) / tr; @@ -473,7 +473,7 @@ int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], dire /**********************************************************************/ void *fields::do_get_array_slice(const volume &where, std::vector components, field_function fun, field_rfunction rfun, void *fun_data, - void *vslice, double omega) { + void *vslice, double frequency) { am_now_working_on(FieldOutput); /***************************************************************/ @@ -510,7 +510,7 @@ void *fields::do_get_array_slice(const volume &where, std::vector com data.rfun = rfun; data.fun_data = fun_data; data.components = components; - data.omega = omega; + data.frequency = frequency; int num_components = components.size(); data.cS = new component[num_components]; data.ph = new cdouble[num_components]; @@ -568,38 +568,38 @@ void *fields::do_get_array_slice(const volume &where, std::vector com /* entry points to get_array_slice */ /***************************************************************/ double *fields::get_array_slice(const volume &where, std::vector components, - field_rfunction rfun, void *fun_data, double *slice, double omega) { - return (double *)do_get_array_slice(where, components, 0, rfun, fun_data, (void *)slice, omega); + field_rfunction rfun, void *fun_data, double *slice, double frequency) { + return (double *)do_get_array_slice(where, components, 0, rfun, fun_data, (void *)slice, frequency); } cdouble *fields::get_complex_array_slice(const volume &where, std::vector components, field_function fun, void *fun_data, cdouble *slice, - double omega) { - return (cdouble *)do_get_array_slice(where, components, fun, 0, fun_data, (void *)slice, omega); + double frequency) { + return (cdouble *)do_get_array_slice(where, components, fun, 0, fun_data, (void *)slice, frequency); } -double *fields::get_array_slice(const volume &where, component c, double *slice, double omega) { +double *fields::get_array_slice(const volume &where, component c, double *slice, double frequency) { std::vector components(1); components[0] = c; return (double *)do_get_array_slice(where, components, 0, default_field_rfunc, 0, (void *)slice, - omega); + frequency); } double *fields::get_array_slice(const volume &where, derived_component c, double *slice, - double omega) { + double frequency) { int nfields; component carray[12]; field_rfunction rfun = derived_component_func(c, gv, nfields, carray); std::vector cs(carray, carray + nfields); - return (double *)do_get_array_slice(where, cs, 0, rfun, &nfields, (void *)slice, omega); + return (double *)do_get_array_slice(where, cs, 0, rfun, &nfields, (void *)slice, frequency); } cdouble *fields::get_complex_array_slice(const volume &where, component c, cdouble *slice, - double omega) { + double frequency) { std::vector components(1); components[0] = c; return (cdouble *)do_get_array_slice(where, components, default_field_func, 0, 0, (void *)slice, - omega); + frequency); } cdouble *fields::get_source_slice(const volume &where, component source_slice_component, diff --git a/src/h5fields.cpp b/src/h5fields.cpp index 652353bdf..1a8cce205 100644 --- a/src/h5fields.cpp +++ b/src/h5fields.cpp @@ -50,7 +50,7 @@ typedef struct { complex *ph; complex *fields; ptrdiff_t *offsets; - double omega; + double frequency; int ninveps; component inveps_cs[3]; direction inveps_ds[3]; @@ -144,7 +144,7 @@ static void h5_output_chunkloop(fields_chunk *fc, int ichnk, component cgrid, iv ptrdiff_t *off = data->offsets; component *cS = data->cS; complex *fields = data->fields, *ph = data->ph; - double omega = data->omega; + double frequency = data->frequency; const component *iecs = data->inveps_cs; const direction *ieds = data->inveps_ds; ptrdiff_t ieos[6]; @@ -175,11 +175,11 @@ static void h5_output_chunkloop(fields_chunk *fc, int ichnk, component cgrid, iv if (cS[i] == Dielectric) { complex tr(0.0,0.0); for (int k = 0; k < data->ninveps; ++k) { - tr += (fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx, omega) + - fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k], omega) + - fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[1 + 2 * k], omega) + + tr += (fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx, frequency) + + fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k], frequency) + + fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[1 + 2 * k], frequency) + fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k] + ieos[1 + 2 * k], - omega)); + frequency)); if (tr == 0.0) tr += 4.0; // default inveps == 1 } fields[i] = (4.0 * data->ninveps) / tr; @@ -187,11 +187,11 @@ static void h5_output_chunkloop(fields_chunk *fc, int ichnk, component cgrid, iv else if (cS[i] == Permeability) { complex tr(0.0,0.0); for (int k = 0; k < data->ninvmu; ++k) { - tr += (fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx, omega) + - fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k], omega) + - fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[1 + 2 * k], omega) + + tr += (fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx, frequency) + + fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k], frequency) + + fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[1 + 2 * k], frequency) + fc->s->get_chi1inv_at_pt(imcs[k], imds[k], idx + imos[2 * k] + imos[1 + 2 * k], - omega)); + frequency)); if (tr == 0.0) tr += 4.0; // default invmu == 1 } fields[i] = (4.0 * data->ninvmu) / tr; @@ -224,7 +224,7 @@ static void h5_output_chunkloop(fields_chunk *fc, int ichnk, component cgrid, iv void fields::output_hdf5(h5file *file, const char *dataname, int num_fields, const component *components, field_function fun, void *fun_data_, int reim, const volume &where, bool append_data, bool single_precision, - double omega) { + double frequency) { am_now_working_on(FieldOutput); h5_output_data data; @@ -270,7 +270,7 @@ void fields::output_hdf5(h5file *file, const char *dataname, int num_fields, data.fun_data_ = fun_data_; /* compute inverse-epsilon directions for computing Dielectric fields */ - data.omega = omega; + data.frequency = frequency; data.ninveps = 0; bool needs_dielectric = false; for (int i = 0; i < num_fields; ++i) @@ -320,23 +320,23 @@ void fields::output_hdf5(h5file *file, const char *dataname, int num_fields, void fields::output_hdf5(const char *dataname, int num_fields, const component *components, field_function fun, void *fun_data_, const volume &where, h5file *file, bool append_data, bool single_precision, const char *prefix, - bool real_part_only, double omega) { + bool real_part_only, double frequency) { bool delete_file; if ((delete_file = !file)) file = open_h5file(dataname, h5file::WRITE, prefix, true); if (real_part_only) { output_hdf5(file, dataname, num_fields, components, fun, fun_data_, 0, where, append_data, - single_precision, omega); + single_precision, frequency); } else { int len = strlen(dataname) + 5; char *dataname2 = new char[len]; snprintf(dataname2, len, "%s%s", dataname, ".r"); output_hdf5(file, dataname2, num_fields, components, fun, fun_data_, 0, where, append_data, - single_precision, omega); + single_precision, frequency); snprintf(dataname2, len, "%s%s", dataname, ".i"); output_hdf5(file, dataname2, num_fields, components, fun, fun_data_, 1, where, append_data, - single_precision, omega); + single_precision, frequency); delete[] dataname2; } if (delete_file) delete file; @@ -357,7 +357,7 @@ static complex rintegrand_fun(const complex *fields, const vec & void fields::output_hdf5(const char *dataname, int num_fields, const component *components, field_rfunction fun, void *fun_data_, const volume &where, h5file *file, bool append_data, bool single_precision, const char *prefix, - double omega) { + double frequency) { bool delete_file; if ((delete_file = !file)) file = open_h5file(dataname, h5file::WRITE, prefix, true); @@ -365,7 +365,7 @@ void fields::output_hdf5(const char *dataname, int num_fields, const component * data.fun = fun; data.fun_data_ = fun_data_; output_hdf5(file, dataname, num_fields, components, rintegrand_fun, (void *)&data, 0, where, - append_data, single_precision, omega); + append_data, single_precision, frequency); if (delete_file) delete file; } @@ -379,27 +379,27 @@ static complex component_fun(const complex *fields, const vec &l } void fields::output_hdf5(component c, const volume &where, h5file *file, bool append_data, - bool single_precision, const char *prefix, double omega) { + bool single_precision, const char *prefix, double frequency) { if (is_derived(int(c))) { - output_hdf5(derived_component(c), where, file, append_data, single_precision, prefix, omega); + output_hdf5(derived_component(c), where, file, append_data, single_precision, prefix, frequency); return; } if (coordinate_mismatch(gv.dim, c)) return; char dataname[256]; - bool has_imag = omega != 0 || (!is_real && c != Dielectric && c != Permeability); + bool has_imag = frequency != 0 || (!is_real && c != Dielectric && c != Permeability); bool delete_file; if ((delete_file = !file)) file = open_h5file(component_name(c), h5file::WRITE, prefix, true); snprintf(dataname, 256, "%s%s", component_name(c), has_imag ? ".r" : ""); output_hdf5(file, dataname, 1, &c, component_fun, 0, 0, where, append_data, single_precision, - omega); + frequency); if (has_imag) { snprintf(dataname, 256, "%s.i", component_name(c)); output_hdf5(file, dataname, 1, &c, component_fun, 0, 1, where, append_data, single_precision, - omega); + frequency); } if (delete_file) delete file; @@ -408,9 +408,9 @@ void fields::output_hdf5(component c, const volume &where, h5file *file, bool ap /***************************************************************************/ void fields::output_hdf5(derived_component c, const volume &where, h5file *file, bool append_data, - bool single_precision, const char *prefix, double omega) { + bool single_precision, const char *prefix, double frequency) { if (!is_derived(int(c))) { - output_hdf5(component(c), where, file, append_data, single_precision, prefix, omega); + output_hdf5(component(c), where, file, append_data, single_precision, prefix, frequency); return; } @@ -421,7 +421,7 @@ void fields::output_hdf5(derived_component c, const volume &where, h5file *file, field_rfunction fun = derived_component_func(c, gv, nfields, cs); output_hdf5(component_name(c), nfields, cs, fun, &nfields, where, file, append_data, - single_precision, prefix, omega); + single_precision, prefix, frequency); } /***************************************************************************/ diff --git a/src/meep.hpp b/src/meep.hpp index f9a9f162e..9ffce676b 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -609,10 +609,10 @@ class structure_chunk { void remove_susceptibilities(); // monitor.cpp - std::complex get_chi1inv_at_pt(component, direction, int idx, double omega = 0) const; - std::complex get_chi1inv(component, direction, const ivec &iloc, double omega = 0) const; - std::complex get_inveps(component c, direction d, const ivec &iloc, double omega = 0) const { - return get_chi1inv(c, d, iloc, omega); + std::complex get_chi1inv_at_pt(component, direction, int idx, double frequency = 0) const; + std::complex get_chi1inv(component, direction, const ivec &iloc, double frequency = 0) const; + std::complex get_inveps(component c, direction d, const ivec &iloc, double frequency = 0) const { + return get_chi1inv(c, d, iloc, frequency); } double max_eps() const; @@ -774,18 +774,18 @@ class structure { void load_chunk_layout(const std::vector &gvs, boundary_region &br); // monitor.cpp - std::complex get_chi1inv(component, direction, const ivec &origloc, double omega = 0, + std::complex get_chi1inv(component, direction, const ivec &origloc, double frequency = 0, bool parallel = true) const; - std::complex get_chi1inv(component, direction, const vec &loc, double omega = 0, + std::complex get_chi1inv(component, direction, const vec &loc, double frequency = 0, bool parallel = true) const; - std::complex get_inveps(component c, direction d, const ivec &origloc, double omega = 0) const { - return get_chi1inv(c, d, origloc, omega); + std::complex get_inveps(component c, direction d, const ivec &origloc, double frequency = 0) const { + return get_chi1inv(c, d, origloc, frequency); } - std::complex get_inveps(component c, direction d, const vec &loc, double omega = 0) const { - return get_chi1inv(c, d, loc, omega); + std::complex get_inveps(component c, direction d, const vec &loc, double frequency = 0) const { + return get_chi1inv(c, d, loc, frequency); } - std::complex get_eps(const vec &loc, double omega = 0) const; - std::complex get_mu(const vec &loc, double omega = 0) const; + std::complex get_eps(const vec &loc, double frequency = 0) const; + std::complex get_mu(const vec &loc, double frequency = 0) const; double max_eps() const; double estimated_cost(int process = my_rank()); @@ -1363,7 +1363,7 @@ class fields_chunk { // monitor.cpp std::complex get_field(component, const ivec &) const; - std::complex get_chi1inv(component, direction, const ivec &iloc, double omega = 0) const; + std::complex get_chi1inv(component, direction, const ivec &iloc, double frequency = 0) const; void backup_component(component c); void average_with_backup(component c); @@ -1551,23 +1551,23 @@ class fields { // low-level function: void output_hdf5(h5file *file, const char *dataname, int num_fields, const component *components, field_function fun, void *fun_data_, int reim, const volume &where, - bool append_data = false, bool single_precision = false, double omega = 0); + bool append_data = false, bool single_precision = false, double frequency = 0); // higher-level functions void output_hdf5(const char *dataname, // OUTPUT COMPLEX-VALUED FUNCTION int num_fields, const component *components, field_function fun, void *fun_data_, const volume &where, h5file *file = 0, bool append_data = false, bool single_precision = false, const char *prefix = 0, - bool real_part_only = false, double omega = 0); + bool real_part_only = false, double frequency = 0); void output_hdf5(const char *dataname, // OUTPUT REAL-VALUED FUNCTION int num_fields, const component *components, field_rfunction fun, void *fun_data_, const volume &where, h5file *file = 0, bool append_data = false, bool single_precision = false, const char *prefix = 0, double = 0); void output_hdf5(component c, // OUTPUT FIELD COMPONENT (or Dielectric) const volume &where, h5file *file = 0, bool append_data = false, - bool single_precision = false, const char *prefix = 0, double omega = 0); + bool single_precision = false, const char *prefix = 0, double frequency = 0); void output_hdf5(derived_component c, // OUTPUT DERIVED FIELD COMPONENT const volume &where, h5file *file = 0, bool append_data = false, - bool single_precision = false, const char *prefix = 0, double omega = 0); + bool single_precision = false, const char *prefix = 0, double frequency = 0); h5file *open_h5file(const char *name, h5file::access_mode mode = h5file::WRITE, const char *prefix = NULL, bool timestamp = false); const char *h5file_name(const char *name, const char *prefix = NULL, bool timestamp = false); @@ -1609,21 +1609,21 @@ class fields { // must eventually be caller-deallocated via delete[]. double *get_array_slice(const volume &where, std::vector components, field_rfunction rfun, void *fun_data, double *slice = 0, - double omega = 0); + double frequency = 0); std::complex *get_complex_array_slice(const volume &where, std::vector components, field_function fun, void *fun_data, - std::complex *slice = 0, double omega = 0); + std::complex *slice = 0, double frequency = 0); // alternative entry points for when you have no field // function, i.e. you want just a single component or // derived component.) - double *get_array_slice(const volume &where, component c, double *slice = 0, double omega = 0); + double *get_array_slice(const volume &where, component c, double *slice = 0, double frequency = 0); double *get_array_slice(const volume &where, derived_component c, double *slice = 0, - double omega = 0); + double frequency = 0); std::complex *get_complex_array_slice(const volume &where, component c, - std::complex *slice = 0, double omega = 0); + std::complex *slice = 0, double frequency = 0); // like get_array_slice, but for *sources* instead of fields std::complex *get_source_slice(const volume &where, component source_slice_component, @@ -1632,7 +1632,7 @@ class fields { // master routine for all above entry points void *do_get_array_slice(const volume &where, std::vector components, field_function fun, field_rfunction rfun, void *fun_data, void *vslice, - double omega = 0); + double frequency = 0); /* fetch and return coordinates and integration weights of grid points covered by an array slice, */ @@ -1679,7 +1679,7 @@ class fields { // that can be passed to eigenmode_amplitude() to get // values of field components at arbitrary points in space. // call destroy_eigenmode_data() to deallocate it when finished. - void *get_eigenmode(double omega_src, direction d, const volume where, const volume eig_vol, + void *get_eigenmode(double frequency, direction d, const volume where, const volume eig_vol, int band_num, const vec &kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, double *kdom = 0, void **user_mdata = 0); @@ -1897,13 +1897,13 @@ class fields { dft_near2far add_dft_near2far(const volume_list *where, const double *freq, size_t Nfreq, int Nperiods = 1); // monitor.cpp - std::complex get_chi1inv(component, direction, const vec &loc, double omega = 0, + std::complex get_chi1inv(component, direction, const vec &loc, double frequency = 0, bool parallel = true) const; - std::complex get_inveps(component c, direction d, const vec &loc, double omega = 0) const { - return get_chi1inv(c, d, loc, omega); + std::complex get_inveps(component c, direction d, const vec &loc, double frequency = 0) const { + return get_chi1inv(c, d, loc, frequency); } - std::complex get_eps(const vec &loc, double omega = 0) const; - std::complex get_mu(const vec &loc, double omega = 0) const; + std::complex get_eps(const vec &loc, double frequency = 0) const; + std::complex get_mu(const vec &loc, double frequency = 0) const; void get_point(monitor_point *p, const vec &) const; monitor_point *get_new_point(const vec &, monitor_point *p = NULL) const; @@ -1983,7 +1983,7 @@ class fields { public: // monitor.cpp std::complex get_field(component c, const ivec &iloc, bool parallel = true) const; - std::complex get_chi1inv(component, direction, const ivec &iloc, double omega = 0, + std::complex get_chi1inv(component, direction, const ivec &iloc, double frequency = 0, bool parallel = true) const; // boundaries.cpp bool locate_component_point(component *, ivec *, std::complex *) const; diff --git a/src/monitor.cpp b/src/monitor.cpp index bd28031fb..3a03638d7 100644 --- a/src/monitor.cpp +++ b/src/monitor.cpp @@ -163,7 +163,7 @@ complex fields_chunk::get_field(component c, const ivec &iloc) const { return 0.0; } -complex fields::get_chi1inv(component c, direction d, const ivec &origloc, double omega, +complex fields::get_chi1inv(component c, direction d, const ivec &origloc, double frequency, bool parallel) const { ivec iloc = origloc; complex aaack = 1.0; @@ -173,53 +173,53 @@ complex fields::get_chi1inv(component c, direction d, const ivec &origlo if (chunks[i]->gv.owns(S.transform(iloc, sn))) { signed_direction ds = S.transform(d, sn); complex val = - chunks[i]->get_chi1inv(S.transform(c, sn), ds.d, S.transform(iloc, sn), omega) * + chunks[i]->get_chi1inv(S.transform(c, sn), ds.d, S.transform(iloc, sn), frequency) * complex(ds.flipped ^ S.transform(component_direction(c), sn).flipped ? -1 : 1,0); return parallel ? sum_to_all(val) : val; } return d == component_direction(c) ? 1.0 : 0; // default to vacuum outside computational cell } -complex fields_chunk::get_chi1inv(component c, direction d, const ivec &iloc, double omega) const { - return s->get_chi1inv(c, d, iloc, omega); +complex fields_chunk::get_chi1inv(component c, direction d, const ivec &iloc, double frequency) const { + return s->get_chi1inv(c, d, iloc, frequency); } -complex fields::get_chi1inv(component c, direction d, const vec &loc, double omega, +complex fields::get_chi1inv(component c, direction d, const vec &loc, double frequency, bool parallel) const { ivec ilocs[8]; double w[8]; complex res(0.0,0.0); gv.interpolate(c, loc, ilocs, w); for (int argh = 0; argh < 8 && w[argh] != 0; argh++) - res += w[argh] * get_chi1inv(c, d, ilocs[argh], omega, false); + res += w[argh] * get_chi1inv(c, d, ilocs[argh], frequency, false); return parallel ? sum_to_all(res) : res; } -complex fields::get_eps(const vec &loc, double omega) const { +complex fields::get_eps(const vec &loc, double frequency) const { complex tr(0.0,0.0); int nc = 0; FOR_ELECTRIC_COMPONENTS(c) { if (gv.has_field(c)) { - tr += get_chi1inv(c, component_direction(c), loc, omega, false); + tr += get_chi1inv(c, component_direction(c), loc, frequency, false); ++nc; } } return complex(nc,0) / sum_to_all(tr); } -complex fields::get_mu(const vec &loc, double omega) const { +complex fields::get_mu(const vec &loc, double frequency) const { complex tr(0.0,0.0); int nc = 0; FOR_MAGNETIC_COMPONENTS(c) { if (gv.has_field(c)) { - tr += get_chi1inv(c, component_direction(c), loc, omega, false); + tr += get_chi1inv(c, component_direction(c), loc, frequency, false); ++nc; } } return complex(nc,0) / sum_to_all(tr); } -complex structure::get_chi1inv(component c, direction d, const ivec &origloc, double omega, +complex structure::get_chi1inv(component c, direction d, const ivec &origloc, double frequency, bool parallel) const { ivec iloc = origloc; for (int sn = 0; sn < S.multiplicity(); sn++) @@ -227,7 +227,7 @@ complex structure::get_chi1inv(component c, direction d, const ivec &ori if (chunks[i]->gv.owns(S.transform(iloc, sn))) { signed_direction ds = S.transform(d, sn); complex val = - chunks[i]->get_chi1inv(S.transform(c, sn), ds.d, S.transform(iloc, sn), omega) * + chunks[i]->get_chi1inv(S.transform(c, sn), ds.d, S.transform(iloc, sn), frequency) * complex((ds.flipped ^ S.transform(component_direction(c), sn).flipped ? -1 : 1),0); return parallel ? sum_to_all(val) : val; } @@ -255,10 +255,10 @@ void matrix_invert(std::complex (&Vinv)[9], std::complex (&V)[9] Vinv[2 + 3 * 2] = 1.0 / det * (V[0 + 3 * 0] * V[1 + 3 * 1] - V[0 + 3 * 1] * V[1 + 3 * 0]); } -complex structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, double omega) const { +complex structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, double frequency) const { complex res(0.0,0.0); if (is_mine()) { - if (omega == 0) + if (frequency == 0) return chi1inv[c][d] ? chi1inv[c][d][idx] : (d == component_direction(c) ? 1.0 : 0); // ----------------------------------------------------------------- // // ---- Step 1: Get instantaneous chi1 tensor ---------------------- @@ -326,7 +326,7 @@ complex structure_chunk::get_chi1inv_at_pt(component c, direction d, int while (my_sus) { if (my_sus->sigma[cc][dd]) { double sigma = my_sus->sigma[cc][dd][idx]; - eps += my_sus->chi1(omega, sigma); + eps += my_sus->chi1(frequency, sigma); } my_sus = my_sus->next; } @@ -334,7 +334,7 @@ complex structure_chunk::get_chi1inv_at_pt(component c, direction d, int // Account for conductivity term if (conductivity[cc][dd]) { double conductivityCur = conductivity[cc][dd][idx]; - eps = std::complex(1.0, (conductivityCur / omega)) * eps; + eps = std::complex(1.0, (conductivityCur / frequency)) * eps; } chi1_tensor[com_it + 3 * dir_int] = eps; } @@ -351,39 +351,39 @@ complex structure_chunk::get_chi1inv_at_pt(component c, direction d, int } complex structure_chunk::get_chi1inv(component c, direction d, const ivec &iloc, - double omega) const { - return get_chi1inv_at_pt(c, d, gv.index(c, iloc), omega); + double frequency) const { + return get_chi1inv_at_pt(c, d, gv.index(c, iloc), frequency); } -complex structure::get_chi1inv(component c, direction d, const vec &loc, double omega, +complex structure::get_chi1inv(component c, direction d, const vec &loc, double frequency, bool parallel) const { ivec ilocs[8]; double w[8]; complex res(0.0,0.0); gv.interpolate(c, loc, ilocs, w); for (int argh = 0; argh < 8 && w[argh] != 0; argh++) - res += w[argh] * get_chi1inv(c, d, ilocs[argh], omega, false); + res += w[argh] * get_chi1inv(c, d, ilocs[argh], frequency, false); return parallel ? sum_to_all(res) : res; } -complex structure::get_eps(const vec &loc, double omega) const { +complex structure::get_eps(const vec &loc, double frequency) const { complex tr(0.0,0.0); int nc = 0; FOR_ELECTRIC_COMPONENTS(c) { if (gv.has_field(c)) { - tr += get_chi1inv(c, component_direction(c), loc, omega, false); + tr += get_chi1inv(c, component_direction(c), loc, frequency, false); ++nc; } } return complex(nc,0) / sum_to_all(tr); } -complex structure::get_mu(const vec &loc, double omega) const { +complex structure::get_mu(const vec &loc, double frequency) const { complex tr(0.0,0.0); int nc = 0; FOR_MAGNETIC_COMPONENTS(c) { if (gv.has_field(c)) { - tr += get_chi1inv(c, component_direction(c), loc, omega, false); + tr += get_chi1inv(c, component_direction(c), loc, frequency, false); ++nc; } } diff --git a/src/mpb.cpp b/src/mpb.cpp index 74ee475ed..cd507e255 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -38,7 +38,7 @@ namespace meep { typedef struct { const double *s, *o; - double omega; + double frequency; ndim dim; const fields *f; } meep_mpb_eps_data; @@ -48,19 +48,19 @@ static void meep_mpb_eps(symmetric_matrix *eps, symmetric_matrix *eps_inv, const meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_; const double *s = eps_data->s; const double *o = eps_data->o; - double omega = eps_data->omega; + double frequency = eps_data->frequency; vec p(eps_data->dim == D3 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1], o[2] + r[2] * s[2]) : (eps_data->dim == D2 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1]) : /* D1 */ vec(o[2] + r[2] * s[2]))); const fields *f = eps_data->f; - eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega)); - eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega)); - eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega)); + eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, frequency)); + eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, frequency)); + eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, frequency)); - ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Y, p, omega)), 0); - ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Z, p, omega)), 0); - ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, omega)), 0); + ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Y, p, frequency)), 0); + ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Z, p, frequency)), 0); + ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, frequency)), 0); /* master_printf("m11(%g,%g) = %g\n", p.x(), p.y(), eps_inv->m00); master_printf("m22(%g,%g) = %g\n", p.x(), p.y(), eps_inv->m11); @@ -102,7 +102,7 @@ typedef struct eigenmode_data { vec center; amplitude_function amp_func; int band_num; - double omega; + double frequency; double group_velocity; } eigenmode_data; @@ -234,7 +234,7 @@ static double dot_product(const mpb_real a[3], const mpb_real b[3]) { } /****************************************************************/ -/* call MPB to get the band_numth eigenmode at freq omega_src. */ +/* call MPB to get the band_numth eigenmode at freq frequency. */ /* */ /* this routine constitutes the first 75% of what was formerly */ /* add_eigenmode_source; it has been split off as a separate */ @@ -251,7 +251,7 @@ static double dot_product(const mpb_real a[3], const mpb_real b[3]) { /* field components at arbitrary points in space. */ /* call destroy_eigenmode_data() to deallocate when finished. */ /****************************************************************/ -void *fields::get_eigenmode(double omega_src, direction d, const volume where, const volume eig_vol, +void *fields::get_eigenmode(double frequency, direction d, const volume where, const volume eig_vol, int band_num, const vec &_kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, double *kdom, void **user_mdata) { @@ -366,7 +366,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c eps_data.o = o; eps_data.dim = gv.dim; eps_data.f = this; - eps_data.omega = omega_src; + eps_data.frequency = frequency; set_maxwell_dielectric(mdata, mesh_size, R, G, meep_mpb_eps, NULL, &eps_data); if (user_mdata) *user_mdata = (void *)mdata; } @@ -400,7 +400,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c // which we automatically pick if kmatch == 0. if (match_frequency && kmatch == 0) { vec cen = eig_vol.center(); - kmatch = omega_src * sqrt(real(get_eps(cen, omega_src)) * real(get_mu(cen, omega_src))); + kmatch = frequency * sqrt(real(get_eps(cen, frequency)) * real(get_mu(cen, frequency))); if (d == NO_DIRECTION) { for (int i = 0; i < 3; ++i) k[i] = dot_product(R[i], kdir) * kmatch; // kdir*kmatch in reciprocal basis @@ -464,7 +464,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c (void *)constraints, W, 3, eigensolver_tol, &num_iters, EIGS_DEFAULT_FLAGS | (am_master() && verbosity > 1 ? EIGS_VERBOSE : 0)); if (verbosity > 0) - master_printf("MPB solved for omega_%d(%g,%g,%g) = %g after %d iters\n", band_num, + master_printf("MPB solved for frequency_%d(%g,%g,%g) = %g after %d iters\n", band_num, G[0][0] * k[0], G[1][1] * k[1], G[2][2] * k[2], sqrt(eigvals[band_num - 1]), num_iters); @@ -486,7 +486,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c if (match_frequency) { // update k via Newton step - double dkmatch = (sqrt(eigvals[band_num - 1]) - omega_src) / vgrp; + double dkmatch = (sqrt(eigvals[band_num - 1]) - frequency) / vgrp; kmatch = kmatch - dkmatch; if (verbosity > 1) master_printf("Newton step: group velocity v=%g, kmatch=%g\n", vgrp, kmatch); @@ -506,7 +506,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c update_maxwell_data_k(mdata, k, G[0], G[1], G[2]); } } while (match_frequency && - fabs(sqrt(eigvals[band_num - 1]) - omega_src) > omega_src * match_tol); + fabs(sqrt(eigvals[band_num - 1]) - frequency) > frequency * match_tol); double eigval = eigvals[band_num - 1]; @@ -529,7 +529,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c if (!am_master()) update_maxwell_data_k(mdata, k, G[0], G[1], G[2]); broadcast(0, (double *)H.data, 2 * H.n * H.p); - if (!match_frequency) omega_src = sqrt(eigval); + if (!match_frequency) frequency = sqrt(eigval); /*--------------------------------------------------------------*/ /*- part 3: do one stage of postprocessing to tabulate H-field */ @@ -578,10 +578,10 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c maxwell_compute_d_from_H(mdata, H, fft_data_E, band_num - 1, 1); - // d_from_H actually computes -omega*D (see mpb/src/maxwell/maxwell_op.c), - // so we need to divide the E-field amplitudes by -omega; we also take this + // d_from_H actually computes -frequency*D (see mpb/src/maxwell/maxwell_op.c), + // so we need to divide the E-field amplitudes by -frequency; we also take this // opportunity to rescale the overall E and H amplitudes to yield unit power flux. - double scale = -1.0 / omega_src, factor = 2.0 / sqrt(vgrp); + double scale = -1.0 / frequency, factor = 2.0 / sqrt(vgrp); cdouble *efield = (cdouble *)fft_data_E, *hfield = (cdouble *)(mdata->fft_data); for (int n = 0; n < NFFT; ++n) { efield[n] *= factor * scale; @@ -610,7 +610,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c edata->center = eig_vol.center(); edata->amp_func = default_amp_func; edata->band_num = band_num; - edata->omega = omega_src; + edata->frequency = frequency; edata->group_velocity = (double)vgrp; if (kdom) { @@ -676,11 +676,11 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d /*--------------------------------------------------------------*/ /* step 1: call MPB to compute the eigenmode */ /*--------------------------------------------------------------*/ - double omega_src = real(src.frequency()); + double frequency = real(src.frequency()); am_now_working_on(MPBTime); global_eigenmode_data = - (eigenmode_data *)get_eigenmode(omega_src, d, where, eig_vol, band_num, kpoint, + (eigenmode_data *)get_eigenmode(frequency, d, where, eig_vol, band_num, kpoint, match_frequency, parity, resolution, eigensolver_tol); finished_working(); @@ -696,7 +696,7 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d global_eigenmode_data->amp_func = A ? A : default_amp_func; src_time *src_mpb = src.clone(); - if (!match_frequency) src_mpb->set_frequency(omega_src); + if (!match_frequency) src_mpb->set_frequency(frequency); /*--------------------------------------------------------------*/ // step 2: add sources whose radiated field reproduces the */ @@ -824,12 +824,12 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in /* dummy versions of class methods for compiling without MPB */ /**************************************************************/ #else // #ifdef HAVE_MPB -void *fields::get_eigenmode(double omega_src, direction d, const volume where, const volume eig_vol, +void *fields::get_eigenmode(double frequency, direction d, const volume where, const volume eig_vol, int band_num, const vec &kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, double *kdom, void **user_mdata) { - (void)omega_src; + (void)frequency; (void)d; (void)where; (void)eig_vol;