diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 8e009dd01..0c51c678c 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -1183,7 +1183,7 @@ Scale the Fourier-transformed fields in `flux` by the complex number `s`. e.g. ` Given a structure, Meep can decompose the Fourier-transformed fields into a superposition of its harmonic modes. For a theoretical background, see [Mode Decomposition](Mode_Decomposition.md). -**`get_eigenmode_coefficients(flux, bands, eig_parity=mp.NO_PARITY, eig_vol=None, eig_resolution=0, eig_tolerance=1e-12, kpoint_func=None, verbose=False)`** +**`get_eigenmode_coefficients(flux, bands, eig_parity=mp.NO_PARITY, eig_vol=None, eig_resolution=0, eig_tolerance=1e-12, kpoint_func=None, verbose=False, direction=mp.AUTOMATIC)`** — Given a flux object and list of band indices, return a `namedtuple` with the following fields: @@ -1194,7 +1194,7 @@ Given a flux object and list of band indices, return a `namedtuple` with the fol The flux object must be created using `add_mode_monitor` (an alias for `add_flux`). `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 the same resolution as Meep, but you can use a higher 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`.) +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). **Note:** for planewaves in homogeneous media, the `kpoints` may *not* necessarily be equivalent to the actual wavevector of the mode. This quantity is given by `kdom`. @@ -1586,7 +1586,7 @@ plt.savefig('sim_domain.png') * `eps_parameters`: a `dict` of optional plotting parameters that override the default parameters for the geometry. - `interpolation='spline36'`: interpolation algorithm used to upsample the pixels. - `cmap='binary'`: the color map of the geometry - - `alpha=1.0`: transparency of geometry + - `alpha=1.0`: transparency of geometry * `boundary_parameters`: a `dict` of optional plotting parameters that override the default parameters for the boundary layers. - `alpha=1.0`: transparency of boundary layers - `facecolor='g'`: color of polygon face diff --git a/python/meep.i b/python/meep.i index 84299e747..a06274d79 100644 --- a/python/meep.i +++ b/python/meep.i @@ -421,7 +421,7 @@ PyObject *_get_dft_array(meep::fields *f, dft_type dft, meep::component c, int n int rank; size_t dims[3]; std::complex *dft_arr = f->get_dft_array(dft, c, num_freq, &rank, dims); - + if (rank==0 || dft_arr==NULL) // this can happen e.g. if component c vanishes by symmetry return PyArray_SimpleNew(0, 0, NPY_CDOUBLE); @@ -492,14 +492,14 @@ kpoint_list get_eigenmode_coefficients_and_kpoints(meep::fields *f, meep::dft_fl int *bands, int num_bands, int parity, double eig_resolution, double eigensolver_tol, std::complex *coeffs, double *vgrp, meep::kpoint_func user_kpoint_func, - void *user_kpoint_data, bool verbose) { + void *user_kpoint_data, bool verbose, meep::direction d) { size_t num_kpoints = num_bands * flux.Nfreq; meep::vec *kpoints = new meep::vec[num_kpoints]; meep::vec *kdom = new meep::vec[num_kpoints]; f->get_eigenmode_coefficients(flux, eig_vol, bands, num_bands, parity, eig_resolution, eigensolver_tol, - coeffs, vgrp, user_kpoint_func, user_kpoint_data, kpoints, kdom, verbose); + coeffs, vgrp, user_kpoint_func, user_kpoint_data, kpoints, kdom, verbose, d); kpoint_list res = {kpoints, num_kpoints, kdom, num_kpoints}; @@ -1353,7 +1353,7 @@ kpoint_list get_eigenmode_coefficients_and_kpoints(meep::fields *f, meep::dft_fl int parity, double eig_resolution, double eigensolver_tol, std::complex *coeffs, double *vgrp, meep::kpoint_func user_kpoint_func, void *user_kpoint_data, - bool verbose); + bool verbose, meep::direction d); PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where, size_t dims[3], bool collapse_empty_dimensions, bool snap_empty_dimensions); diff --git a/python/simulation.py b/python/simulation.py index 0ab46fd08..373fd6bcf 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -148,7 +148,7 @@ class Identity(Symmetry): class Volume(object): def __init__(self, center=Vector3(), size=Vector3(), dims=2, is_cylindrical=False, vertices=[]): - + if len(vertices) == 0: self.center = center self.size = size @@ -164,7 +164,7 @@ def __init__(self, center=Vector3(), size=Vector3(), dims=2, is_cylindrical=Fals z_size = 0 if z_list.size == 1 else np.abs(np.diff(z_list)[0]) self.size = Vector3(x_size,y_size,z_size) - + self.dims = dims v1 = self.center - self.size.scale(0.5) @@ -1934,13 +1934,15 @@ def get_dft_array_metadata(self, dft_cell=None, vol=None, center=None, size=None center=center, size=size, collapse=True) def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, eig_vol=None, - eig_resolution=0, eig_tolerance=1e-12, kpoint_func=None, verbose=False): + eig_resolution=0, eig_tolerance=1e-12, kpoint_func=None, verbose=False, direction=mp.AUTOMATIC): if self.fields is None: raise ValueError("Fields must be initialized before calling get_eigenmode_coefficients") if eig_vol is None: eig_vol = flux.where else: eig_vol = self._volume_from_kwargs(vol=eig_vol) + if direction is None or direction == mp.AUTOMATIC: + direction = flux.normal_direction num_bands = len(bands) coeffs = np.zeros(2 * num_bands * flux.Nfreq, dtype=np.complex128) @@ -1957,7 +1959,8 @@ def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, eig_v coeffs, vgrp, kpoint_func, - verbose + verbose, + direction ) return EigCoeffsResult(np.reshape(coeffs, (num_bands, flux.Nfreq, 2)), vgrp, kpoints, kdom) diff --git a/python/tests/oblique_source.py b/python/tests/oblique_source.py index e53cde5db..a2ff9dd57 100644 --- a/python/tests/oblique_source.py +++ b/python/tests/oblique_source.py @@ -14,13 +14,15 @@ def test_waveguide_flux(self): rot_angles = range(0,60,20) # rotation angle of waveguide, CCW around z-axis fluxes = [] + coeff_fluxes = [] for t in rot_angles: rot_angle = math.radians(t) + kpoint = mp.Vector3(math.cos(rot_angle),math.sin(rot_angle),0) sources = [mp.EigenModeSource(src=mp.GaussianSource(1.0,fwidth=0.1), size=mp.Vector3(y=10), center=mp.Vector3(x=-3), direction=mp.NO_DIRECTION, - eig_kpoint=mp.Vector3(math.cos(rot_angle),math.sin(rot_angle),0), + eig_kpoint=kpoint, eig_band=1, eig_parity=mp.ODD_Z, eig_match_freq=True)] @@ -41,11 +43,19 @@ def test_waveguide_flux(self): sim.run(until_after_sources=100) + res = sim.get_eigenmode_coefficients(tran, [1], + eig_parity=mp.ODD_Z, direction=mp.NO_DIRECTION, + kpoint_func=lambda f,n: kpoint) + fluxes.append(mp.get_fluxes(tran)[0]) + coeff_fluxes.append(abs(res.alpha[0,0,0])**2) print("flux:, {:.2f}, {:.6f}".format(t,fluxes[-1])) + print("coef_flux:, {:.2f}, {:.6f}".format(t,coeff_fluxes[-1])) self.assertAlmostEqual(fluxes[0], fluxes[1], places=0) self.assertAlmostEqual(fluxes[1], fluxes[2], places=0) + for i in range(3): + self.assertAlmostEqual(fluxes[i], coeff_fluxes[i], places=0) # self.assertAlmostEqual(fluxes[0], fluxes[2], places=0) # sadly the above line requires a workaround due to the diff --git a/src/meep.hpp b/src/meep.hpp index 647ad4a8d..d4f51a198 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1669,6 +1669,11 @@ class fields { double eigensolver_tol, std::complex amp, std::complex A(const vec &) = 0); + void get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, int *bands, int num_bands, + int parity, double eig_resolution, double eigensolver_tol, + std::complex *coeffs, double *vgrp, + kpoint_func user_kpoint_func, void *user_kpoint_data, + vec *kpoints, vec *kdom, bool verbose, direction d); void get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, int *bands, int num_bands, int parity, double eig_resolution, double eigensolver_tol, std::complex *coeffs, double *vgrp, diff --git a/src/mpb.cpp b/src/mpb.cpp index 3e7ba92bc..c60fe0dc9 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -57,7 +57,7 @@ static void meep_mpb_eps(symmetric_matrix *eps, symmetric_matrix *eps_inv, const eps_inv->m00 = f->get_chi1inv(Ex, X, p, omega); eps_inv->m11 = f->get_chi1inv(Ey, Y, p, omega); eps_inv->m22 = f->get_chi1inv(Ez, Z, p, omega); - + ASSIGN_ESCALAR(eps_inv->m01, f->get_chi1inv(Ex, Y, p, omega), 0); ASSIGN_ESCALAR(eps_inv->m02, f->get_chi1inv(Ex, Z, p, omega), 0); ASSIGN_ESCALAR(eps_inv->m12, f->get_chi1inv(Ey, Z, p, omega), 0); @@ -746,15 +746,12 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in double eigensolver_tol, std::complex *coeffs, double *vgrp, kpoint_func user_kpoint_func, void *user_kpoint_data, vec *kpoints, vec *kdom_list, - bool verbose) { + bool verbose, direction d) { double freq_min = flux.freq_min; double dfreq = flux.dfreq; int num_freqs = flux.Nfreq; - direction d = flux.normal_direction; bool match_frequency = true; - if (d == NO_DIRECTION) abort("cannot determine normal direction in get_eigenmode_coefficients"); - if (flux.use_symmetry && S.multiplicity() > 1 && parity == 0) abort("flux regions for eigenmode projection with symmetry should be created by " "add_mode_monitor()"); @@ -862,7 +859,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in double eigensolver_tol, std::complex *coeffs, double *vgrp, kpoint_func user_kpoint_func, void *user_kpoint_data, vec *kpoints, vec *kdom, - bool verbose) { + bool verbose, direction d) { (void)flux; (void)eig_vol; (void)bands; @@ -877,6 +874,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in (void)user_kpoint_data; (void)kdom; (void)verbose; + (void) d; abort("Meep must be configured/compiled with MPB for get_eigenmode_coefficient"); } @@ -904,4 +902,15 @@ vec get_k(void *vedata) { #endif // HAVE_MPB +/* compatibility wrapper routine that passes the default flux.normal_direction to the eigensolver + (we pass NO_DIRECTION to use the kpoint direction instead, for oblique sources). */ +void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, int *bands, + int num_bands, int parity, double eig_resolution, + double eigensolver_tol, std::complex *coeffs, + double *vgrp, kpoint_func user_kpoint_func, + void *user_kpoint_data, vec *kpoints, vec *kdom, bool verbose) { + get_eigenmode_coefficients(flux, eig_vol, bands, num_bands, parity, eig_resolution, eigensolver_tol, + coeffs, vgrp, user_kpoint_func, user_kpoint_data, kpoints, kdom, verbose, flux.normal_direction); +} + } // namespace meep